// The Magma code in this file computes the zeta function of the 
// surface Y(mu(3)) over the field F_q, for a prime power q. It is
// used in the file "orbit-counting-A5-example.txt".


function AConst(F,C)
// Given elliptic curve F over a finite field with 
// 2-torsion point C, returns the constant A:=g0*g1.
	if #F eq 2 then
		F2:=BaseExtend(F,2);
		C2:=F2!C;
		return BaseRing(F)!AConst(F2,C2);
	end if;
	repeat 
		f:=Random(F);
	until f ne F!0 and f ne C;
	return (f[1]-C[1])*((f+C)[1]-C[1]);
end function;


function countFiber(F, C, f, a)
// Given elliptic curve F over a finite field with 2-torsion point C
// and arbitrary point f, and a in P^6, counts the number of points in the 
// fiber of Y(a) over f. If the fiber is nonsingular,
// also returns the fiber as a HyperellipticCurve.
	P<Z>:= PolynomialRing(BaseRing(F));
	if f ne F!0 and f ne C then
		g0:= f[1]-C[1];
		h0:= f[2];
		g1:= (f+C)[1]-C[1];
		h1:= (f+C)[2];
		phi0:= g0*Z^6 + g1;
		phi1:= g0^2*Z^6 + g1^2;
		phi2:= g0*Z^5 + g1*Z;
		phi3:= h0*Z^5 + h1*Z;
		phi4:= Z^4 + Z^2;
		phi5:= g0*Z^4 + g1*Z^2;
		phi6:= Z^3;
		q:= a[1]*phi0 + a[2]*phi1 + a[3]*phi2 + a[4]*phi3
			+ a[5]*phi4 + a[6]*phi5 + a[7]*phi6;
		if IsHyperellipticCurveOfGenus(2,[q,0]) then
			H:=HyperellipticCurve(q);
			return #H, H;
		else
			// printf "Singular fiber for f=%o\n", f;
			// print q;
			if Degree(q) eq 6 then
				lead:=BaseRing(F)!Coefficient(q,6);
				// 2 points at infinity if leading coefficient is square; zero if not							
				count:=#AllRoots(lead,2);
			else
				// 1 point at infinity
				count:=1; 
			end if;
			R:=PolynomialRing(BaseRing(F));
			for x in BaseRing(F) do
				count +:= #AllRoots(Evaluate(R!q,x),2);
			end for;
			return count;
		end if;
	else
		A:=AConst(F,C);
		q:=a[1]*A + a[2]*Z^6 + a[4]*(Z^5-A*Z)
			+ a[5]*Z^2 + a[6]*Z^4;
		if IsHyperellipticCurveOfGenus(2,[q,0]) then
			H:=HyperellipticCurve(q);
			return #H, H;
		else
			// printf "Singular fiber for f=%o\n", f;
			// print q;
			if Degree(q) eq 6 then
				lead:=BaseRing(F)!Coefficient(q,6);
				// 2 points at infinity if leading coefficient is square; zero if not
				count:=#AllRoots(lead,2);  
			else
				// 1 point at infinity
				count:=1;
			end if;
			R:=PolynomialRing(BaseRing(F));
			for x in BaseRing(F) do
				count +:= #AllRoots(Evaluate(R!q,x),2);
			end for;
			return count;
		end if;
	end if;
end function;


function countPoints(F,C,a)
// Given elliptic curve F over a finite field with 2-torsion point C and
// a in P^6, counts the number of points of Y(a) defined over the same 
// field as F.
	count:= 0;
	for f in Points(F) do
		count+:= countFiber(F,C,f,a);
	end for;
	printf "Point total is %o\n", count;
	return count;
end function;


function countExtension(F,C,a,n)
// Given elliptic curve F over F_q, 2-torsion point C, and a in P^6, counts 
// the points of Y(a) defined over F_{q^n}.
	Fn:=BaseExtend(F,n);
	Cn:=Fn!C;
	printf "Now counting points over extension of degree %o\n", n;
	return countPoints(Fn,Cn,a);
end function;



procedure computeZeta4(q)
	tau := 3;
	F := EllipticCurve([GF(q)!0,(1-6*tau^2-3*tau^4),0,16*tau^6,0]);
	C := F![0,0,1];
	P6 := ProjectiveSpace(Rationals(),6);
	a := P6![
		-4*tau^2*(1 + tau)^6, 
		(1 + 6*tau + tau^2)^2, 
		4*(-1 + tau)^2*tau^2*(1 + 14*tau + 34*tau^2 + 14*tau^3 + tau^4),
     	0, 
     	4*tau^4*(1 + 12*tau + 40*tau^2 + 164*tau^3 + 334*tau^4 + 164*tau^5 + 40*tau^6 + 12*tau^7 + tau^8), 
     	-8*tau^3*(3 + 28*tau + 34*tau^2 + 28*tau^3 + 3*tau^4), 
     	-8*(-1 + tau)^2*tau^4*(1 + 14*tau + 71*tau^2 + 84*tau^3 + 71*tau^4 + 14*tau^5 + tau^6)
     	];
	print "a-value: ", a;
	time N:= [countExtension(F,C,a,n): n in [1..4]];
	Z:= ZetaFunction(F);
	P1:= Numerator(Z);
	print "L-Polynomial of elliptic curve: ", P1;
	d:= Discriminant(P1);
	Q<t>:= PolynomialRing(QuadraticField(d));
	alpha1:= 1/(Roots(Q!P1)[1][1]);
	alpha2:= 1/(Roots(Q!P1)[2][1]);
	// If the A5 singularities are defined over GF(q), the following quantity will be 1; otherwise they are defined over GF(q^2) and the quantity will be -1.
	A5sign := -1+#AllRoots(GF(q)!((1-tau)*(1+3*tau)),2);
	print "A5sign = ",A5sign;
	print "Factor from fiber, taut. div., and two A5 singularities: ";
	Factorization((1-q*t)^7*(1-A5sign*q*t)^5);
	if Norm(alpha1) ne q or Norm(alpha2) ne q
		then print "Roots have wrong absolute value";
	end if;
	S:= [N[r] + 10*(1+(A5sign)^r)/2*q^r - (1+8*q^r+5*(A5sign*q)^r+q^(2*r)) + (q^r+1)*alpha1^r + (q^r+1)*alpha2^r: r in [1..4]];
	// print S;
	B:=[-S[1]];
	for k:= 2 to 4 do
		Bk:=-(1/k) * ((&+[B[k-j]*S[j]: j in [1..k-1]])+S[k]);
		Append(~B, Bk);
	end for;
	// print B;
	print "Assuming sign of functional equation is 1:";
	sign:= -1;
	Bpos:= B;
	for k:= 5 to 8 do
		Append(~Bpos,sign*q^(2*k-9)*Bpos[9-k]);
	end for;
	Append(~Bpos,sign*q^9);
	P2pos:= (1-q*t)^8*(1-A5sign*q*t)^5*(1 + &+[Bpos[k]*t^k: k in [1..9]]);
	print P2pos;
	Factorization(P2pos);
	print "Assuming sign of functional equation is -1:";
	sign:= 1;
	Bneg:= B;
	for k:= 5 to 8 do
		Append(~Bneg,sign*q^(2*k-9)*Bneg[9-k]);
	end for;
	Append(~Bneg,sign*q^9);
	P2neg:= (1-q*t)^8*(1-A5sign*q*t)^5*(1 + &+[Bneg[k]*t^k: k in [1..9]]);
	print P2neg;
	Factorization(P2neg);
end procedure;


// Use the command: computeZeta4(q);