\\ The goal of this Magma code is to figure out, for a given small prime number p,

\\ if there is a quotient of $J_0(p)^+$ with rank=dimension, and dimension $\geq 2$,

\\ and similarly with $J_0(p^2)^{+,\textrm{new}}$.



\\ This function takes $f$ and $N$ where $f$ is in $S_2(\Gamma_0(N))^+$ is a normalised eigenform (with coefficients 

\\ embedded in $\C$) and tests, using the classical formula, if the derivative $L'(f,1)$ is assured to be nonzero. 



\\ The next two functions compute the genuses of $X_0^+(p)$ and $X_{\rm{ns}}^+(p)$, through well-known formulas.



GenusPlusPartX0 := function(p)

g := 1/12*(p-6-3*LegendreSymbol(-1,p)-4*LegendreSymbol(-3,p));

if (p mod 4 eq 3) then
	V:= ClassNumber(-p) + ClassNumber(-4*p);

else
	V:= ClassNumber(-4*p);
end if;
return ((g+1)/2 - V/4);

end function;



GenusNonSplit := function(p)

g := 1/24*(p^2-10*p+23+6*LegendreSymbol(-1,p)+4*LegendreSymbol(-3,p));

return g;

end function;



\\ This formula tests if the central value of the derivative of the modular form is nonzero (with precision $10^{-3}$, largely sufficient).



TestNonZeroDer := function(f,N)

Re:= RealField(30);

range:=Floor(N^(3/4)/(2*Pi(Re)));

approxL:=&+[Coefficient(f,n)/n*ExponentialIntegralE1(2*Pi(Re)*n/Sqrt(N)): n in [1..range]];

absval:=AbsoluteValue(approxL);

return  (absval ge 1/1000);

end function;


\\ The next two functions compute approximate values for the $L'(f,1)$ in the nonsplit case, and if 2 of them are nonzero, quadratic Chabauty condition is satisfied. 
\\ It gives a complete answer, but is lengthy to compute for large levels, hence the preliminary work in Sage.


NonZeroDerNonSplit := function(p)

NFs:=Newforms(CuspForms(p^2));

dim:=0;

check:=false;

NumberOrbits:=#NFs;

for k in [1..NumberOrbits] do
 
f:=NFs[k][1];
n:=#NFs[k];

g:=ComplexEmbeddings(f)[1,1];

L:=LSeries(g);
	
a:=CheckFunctionalEquation(L);
	
sign:=Sign(L);
	
print sign;

if (AbsoluteValue(sign+1) lt (0.00001)) then
 test := TestNonZeroDer(g,p^2);

if test then
 dim:=dim+n;
 
else print "likely vanishing at", p;

end if;

if (dim ge 2) then
 check:=true;

break;	
end if;
	
end if;

end for;

return check;

end function;





FullNonZeroDerNonSplit := function(p)

NFs:=Newforms(CuspForms(p^2));

dim:=0; 

genus:=GenusNonSplit(p);

check:=false;

NumberOrbits:=#NFs;

for k in [1..NumberOrbits] do
 f:=NFs[k][1];
n:=#NFs[k];

g:=ComplexEmbeddings(f)[1,1];

L:=LSeries(g);
	
a:=CheckFunctionalEquation(L); 
sign:=Sign(L);
	
if (AbsoluteValue(sign+1) lt (0.0000001)) then
	test := TestNonZeroDer(g,p^2);

if test then
 dim:=dim+n;
else print "likely vanishing at", p;

end if;
		
if (dim eq genus) then
	check:=true;

break;

end if;
	
end if;

end for;

return check;

end function;



\\ This function checks that the weaker quadratic Chabauty condition is satisfied in ranges of primes, using the previous function.


FullQuadChabCheckNonSplit := function(M,N)

exceptions:=[];

for p in [M..N] do
 
if IsPrime(p) then
 if FullNonZeroDerNonSplit(p) then print "Whole jacobian is good for p=",p;

else Include(~exceptions,p);

end if;
	
end if;

end for;

return exceptions;

end function;




\\ This function checks that the $L'(f,1)$ are nonzero for a level $p$ prime, the following one does it in a range of primes.

NonZeroDer := function(p)

NFs:=Newforms(CuspForms(p));

dim:=0;

check:=false;

NumberOrbits:=#NFs;

for k in [1..NumberOrbits] do
 f:=NFs[k];
n:=#f;
	
sign:=Coefficient(f[1],p);

if (sign eq (-1)) then
	print "degree of field is", n;

g:=ComplexEmbeddings(f[1])[1,1];
test := TestNonZeroDer(g,p);

if test then	dim:=dim+n;
else print "likely vanishing at", p;

end if;

if (dim ge 2) then check:=true;
break;

end if;
	
end if;

end for;

return check;

end function;



FullNonZeroDer := function(p)

NFs:=Newforms(CuspForms(p));

genus:=GenusPlusPartX0(p);

dim:=0;

NumberOrbits:=#NFs;

for k in [1..NumberOrbits] do
 f:=NFs[k];
n:=#f;
	
sign:=Coefficient(f[1],p);
 
if (sign eq (-1)) then
 g:=ComplexEmbeddings(f[1])[1,1];

test := TestNonZeroDer(g,p);

if test then	dim:=dim+n;

else 
	print "likely vanishing at", p;

end if;
	
end if;

end for;

return (dim eq genus);

end function;




\\ The next auxiliary functions checks whether the considered jacobians are simple or not, and if there is a large chunk apart from elliptic curves.


SimpleJacobian := function(p)

NFs:=Newforms(CuspForms(p));

print "newforn basis computed for p=",p;

dim:=0;

check:=false; 

genus:=GenusPlusPartX0(p);

NumberOrbits:=#NFs;

for k in [1..NumberOrbits] do
 f:=NFs[k];
n:=#f;
	
sign:=Coefficient(f[1],p);

if (sign eq (-1)) then
	check:=(n eq genus);

break;
	
end if;

end for;

return check;

end function;



FactorsJacobianUptoEllCurve:=function(p)

NFs:=Newforms(CuspForms(p));

NumberOrbits:=#NFs;

listdimfactors:=[];

for k in [1..NumberOrbits] do
	f:=NFs[k];
n:=#f;
	
sign:=Coefficient(f[1],p);
if (sign eq (-1)) and n gt 1 then
 Include(~listdimfactors,n);

end if;

end for;

return listdimfactors;

end function;



FactorsJacobianUptoEllCurveRange := function(M,N)

count:=0;

others:=0;

for p in [M..N] do
 if IsPrime(p) then l:=FactorsJacobianUptoEllCurve(p);

if #l eq 1 then
	count:=count+1;

else	others:=others+1;

end if;	
	
end if;

d for;

return count,"simple jacobians outside elliptic curves and",others,"nonsimple ones";

end function;



CheckSimpleJacobian:=function(M,N)

exceptions:=[];

count:=0;

for p in [M..N] do
 if IsPrime(p) then
	if SimpleJacobian(p) then
 count:=count+1;

else	Include(~exceptions,p);

end if;
	
end if;

end for;

print count,"simple jacobians found";
print "exceptions are",exceptions;

end function;





\\ This function computes the genus of X_0(p)^+ given known formulae.



\\ This function checks that quadratic Chabauty hypothesis holds between M and N (when it can, i.e. when the genus is at least 2). All counterexamples are listed during the process.




CheckQuadChabFull := function(M,N)

l:=[];

for p in [M..N] do

if IsPrime(p) then if GenusPlusPartX0(p) ge 2 then if (not FullNonZeroDer(p)) then
 Include(~l,p);
end if;
	
end if;

end if;

end for;

return "list of exceptions in this range: ",l;

end function;




FirstCounterExampleQuadChabFull := function(M,N)

check:=true;

exception:=1;

for p in [M..N] do

if IsPrime(p) then
 if GenusPlusPartX0(p) ge 2 then if (not FullNonZeroDer(p)) then exception:=p;

break;
else print "full jacobian is valid for p=",p;

end if;
	
end if;

end if;

end for;

if check then
	return "no exception";

else
	return "possible exception for p= ",exception;

end if;

end function;



\\ The last functions are made to check the last possible exceptions (see end of section 5 in paper), by computing the contribution of the elliptic curves factors to the moments and checking that it is never enough in those cases.

CheckQuadChabViaEllCurves := function(M,N)

D:=CremonaDatabase();

Re:=RealField(20);

exceptions:=[];

still_exceptions:=[];

for p in [M..N] do if IsPrime(p) then
 n:=NumberOfIsogenyClasses(D,p);

numberellcurves:=0;

for i in [1..n] do
 E:=EllipticCurve(D,p,i,1);

if (RootNumber(E) eq -1) then
	r:=AnalyticRank(E);

if (r eq 1) then
 numberellcurves:=numberellcurves+1;

Ecand:=E;

else 
	print "Elliptic curve with higher rank",r,"and conductor ",p;

end if;
end if;
end for;
if (numberellcurves eq 1) then
	if (TraceOfFrobenius(Ecand,2) eq 1) or ((TraceOfFrobenius(Ecand,2) eq 0) and (p lt 3001)) then 
	r, Lvalue:=AnalyticRank(E);

moddegree:=ModularDegree(E);

omega1:=Periods(E)[1];

omega2:=Periods(E)[2];
volE:=AbsoluteValue(omega1*Imaginary(omega2));
quotient:=4*Pi(Re)^2*Lvalue/(moddegree*volE);

if (quotient le 4/5) then	Include(~exceptions,p);

else	Include(~still_exceptions,p);

end if;
end if;
end if;
	
end if;

end for;

return still_exceptions,exceptions;

end function;




CheckQuadChabViaEllCurves2 := function(M,N)

D:=CremonaDatabase();

Re:=RealField(20);

exceptions:=[];

still_exceptions:=[];

for p in [M..N] do if IsPrime(p) then
	n:=NumberOfIsogenyClasses(D,p);

numberellcurves:=0;
for i in [1..n] do
 E:=EllipticCurve(D,p,i,1);

if (RootNumber(E) eq -1) then
	r:=AnalyticRank(E);

if (r eq 1) then
	numberellcurves:=numberellcurves+1;

Ecand:=E;

else print "Elliptic curve with higher rank",r,"and conductor ",p;

end if;
end if;
end for;

if (numberellcurves eq 1) then
	r, Lvalue:=AnalyticRank(E);

moddegree:=ModularDegree(E);

omega1:=Periods(E)[1];

omega2:=Periods(E)[2];
volE:=AbsoluteValue(omega1*Imaginary(omega2));

quotient:=4*Pi(Re)^2*Lvalue/(moddegree*volE);
if (quotient le 4/5) then	Include(~exceptions,p);

else	Include(~still_exceptions,p);
end if;	
end if;
	
end if;

end for;

return still_exceptions,exceptions;

end function;




CheckQuadChabViaEllCurvesBasic := function(M,N)

D:=CremonaDatabase();

Re:=RealField(20);

exceptions:=[];

for p in [M..N] do	if IsPrime(p) then	n:=NumberOfIsogenyClasses(D,p);

numberellcurves:=0;
for i in [1..n] do
	E:=EllipticCurve(D,p,i,1);

if (RootNumber(E) eq -1) then
	r:=AnalyticRank(E);

Ecand:=E;

if (r eq 1) then
	numberellcurves:=numberellcurves+1;

Ecand:=E;

else 
	print "Elliptic curve with higher rank",r,"and conductor ",p;

end if;
end if;
end for;
if (numberellcurves eq 1) then
	if (TraceOfFrobenius(Ecand,2) eq 1) or ((TraceOfFrobenius(Ecand,2) eq 0) and (p lt 3001)) then 
	Include(~exceptions,p);

end if;

end if;
end if;

end for;

return exceptions;

end function;





CheckQuadChabNonSplitViaEllCurves := function(M,N)

D:=CremonaDatabase();

Re:=RealField(20);

exceptions:=[];

still_exceptions:=[];

for p in [M..N] do
	if IsPrime(p) then
	n:=NumberOfIsogenyClasses(D,p^2);

numberellcurves:=0;

for i in [1..n] do	E:=EllipticCurve(D,p^2,i,1);

if (RootNumber(E) eq -1) then	r:=AnalyticRank(E);

if (r eq 1) then	numberellcurves:=numberellcurves+1;

Ecand:=E;

else 	print "Elliptic curve with higher rank",r,"and conductor ",p^2;

end if;
end if;
end for;

if (numberellcurves eq 1) then
	if (TraceOfFrobenius(Ecand,2) eq 1) or ((TraceOfFrobenius(Ecand,2) eq 0) and (p lt 71))  then r, Lvalue:=AnalyticRank(Ecand);

moddegree:=ModularDegree(Ecand);

omega1:=Periods(Ecand)[1];
	
omega2:=Periods(Ecand)[2];
volE:=AbsoluteValue(omega1*Imaginary(omega2));

quotient:=4*Pi(Re)^2*Lvalue/(moddegree*volE);
if (quotient le 1) then		Include(~exceptions,p);
else	Include(~still_exceptions,p);

end if;
end if;
end if;
end if;

end for;

return still_exceptions,exceptions;

end function;





CheckQuadChabNonSplitViaEllCurvesBasic := function(M,N)

D:=CremonaDatabase();

Re:=RealField(20);

exceptions:=[];

for p in [M..N] do	if IsPrime(p) then
	n:=NumberOfIsogenyClasses(D,p^2);

numberellcurves:=0;

for i in [1..n] do
	E:=EllipticCurve(D,p^2,i,1);

if (RootNumber(E) eq -1) then
	r:=AnalyticRank(E);

if (r eq 1) then
	numberellcurves:=numberellcurves+1;
Ecand:=E;
else	print "Elliptic curve with higher rank",r,"and conductor ",p^2;

end if;
end if;
end for;
if (numberellcurves eq 1) then	if (TraceOfFrobenius(Ecand,2) eq 1)  or ((TraceOfFrobenius(Ecand,2) eq 0) and (p lt 71)) then 
	Include(~exceptions,p);
end if;
end if;
	
end if;

end for;

return exceptions;

end function;



CheckQuadChabPrime := function(M,N)

exceptions:=[];

for p in [M..N] do
	if IsPrime(p) then
	if GenusPlusPartX0(p) ge 2 then	if (not NonZeroDer(p)) then 
exceptions:=exceptions+[p];

print "exception found for p=", p;

else	print "QCH holds for ",p;
end if;
end if;
	
end if;

end for;

return exceptions;

end function;



CheckQuadChabPrimeList := function(l)

exceptions:=[];

for p in l do
	if (not NonZeroDer(p)) then	exceptions:=exceptions+[p];

print "exception found for p=", p;

else	print "QCH holds for ",p;

end if;

end for;

return exceptions;

end function;