
// Linear combination needed in Euclidean algorithm routine.
// It returns a,b such that a*f+b*g has lower degree than g.

function Red(f,g)

	x := Parent(f).1;	
	L := [LeadingCoefficient(g,1), (-x^(Degree(f,1)-Degree(g,1))*LeadingCoefficient(f,1))];
	F := L[1]*f + L[2]*g;
	while Degree(F,1) ge Degree(g,1) do;
		L := [L[1]*LeadingCoefficient(g,1),L[2]*LeadingCoefficient(g,1)-x^(Degree(F,1)-Degree(g,1))*LeadingCoefficient(F,1)];
		F := LeadingCoefficient(g,1)*F + (-x^(Degree(F,1)-Degree(g,1))*LeadingCoefficient(F,1))*g;
	end while;
	
	return L[1],L[2],g,F;	

end function;

// Euclidean algorithm in the first variable x, return p,q such that p*f+q*g = 1.
// We assume that degree(f) >= degree(g).

function EuclidX(f,g)

	F := f;
	G := g;
	L := [];
	
	while Degree(G,1) gt 0 do
		a,b,F,G := Red(F,G);
		Append(~L,[a,b]);
	end while;
	
	l := #L;
	if l eq 0 then // Don't know if this ever happens when you come from main function, but oh well.
		A := 0;
		B := 1;
	else
		A := L[l][1];
		B := L[l][2];
	end if;
	
	for i := l-1 to 1 by -1 do
		a := A; 
		b := B; 
		A:= b*L[i][1];
		B := b*L[i][2]+a;
	end for;
	
	R<x,t> := FieldOfFractions(Parent(f));
		
	assert (R!A/G)*f+(R!B/G)*g eq 1;
	return R!A/G,R!B/G;

end function;

// Reduction algorithm in DeRham cohomology. Algorithm subtracts exact forms until we are left with a linear combination of omega and xomega.

function KedlayaRed(f,P)

	x := Parent(f).1;
	deg := Degree(P,1);
	while Degree(Numerator(f),1) ge deg-1 do
		D := Degree(Numerator(f),1);
		C := LeadingCoefficient(Numerator(f),1)/Denominator(f);
		F := C/(D-(deg-1)+deg/2)*x^(D-(deg-1));
		f := f - (Derivative(P,1)*F/2 + Derivative(F,1)*P);
	end while;
	
	return f;

end function;

// This function returns the Gauss-Manin connection matrix on the family of elliptic curves defined by P (assumed monic).

function GaussManinMatrix(P)

	assert LeadingCoefficient(P) eq 1;
	d := Degree(P,1);
	g := Floor((d-1)/2);
	// Find the genus of the hyperelliptic curve.
	genus := Floor((d-1)/2);

	R<x,t> := PolynomialRing(Rationals(),2);
	P := R!P;
	Px := Derivative(P,1);
	Pt := Derivative(P,2);
		
	// Define A,B as appearing in Kedlaya's AWS notes on p.23
	A,B := EuclidX(P,Px);
	Ax := Derivative(A,1);
	At := Derivative(A,2);
	Bx := Derivative(B,1);
	Bt := Derivative(B,2);
	
	// First calculate the vector [D x^i omega]_i.
	L := [Bx*Pt - A*Pt/2 - At*P - Bt*Px];
	for i := 1 to d-2 do
		Append(~L,x^i*L[1]+i*x^(i-1)*B*Pt);
	end for;
	
	// Now reduce it by subtracting exact forms, in basis x^idx/y
	D := [];
	for i := 1 to d-1 do
		Append(~D,KedlayaRed(L[i],P));
	end for;

	// Calculate the matrix of GM wrt basis x^idx/y.
	M := ZeroMatrix(FieldOfFractions(R),d-1,d-1);
	for i := 1 to d-1 do
		for j := 1 to d-1 do
			M[j,i] := (Coefficient(Numerator(D[i]),1,j-1)/Denominator(D[i]));
		end for;
	end for;
	
	// Calculate the matrix of GM wrt some to be determined basis (residue at infinity = 0) deduced from x^idx/y.
	if IsEven(d) then
		// First compute the residues of x^idx/y at infinity for d even.
		S<W> := LaurentSeriesRing(R);
		AW := S!-2;
		for i:= 1 to 40 do
			C := 0;
			if d-i ge 0 then
				C := Coefficient(P,1,d-i);
			end if;
			AW := AW - (Coefficient(-AW^2,i)+C)*W^(i)/2;
		end for;
		
		// Now calculate the residues of x^idx/y at (infty+).
		ResList:=[];
		Series := -1/(AW+1);
		for i := 1 to d-1 do
			Append(~ResList,Coefficient(Series,i-g-1));
		end for;
		// Compute smallest index with nonzero residue.
		ind := 1;
		while ResList[ind] eq 0 do
			ind := ind + 1;
		end while;
		
		// Calculate the vector [D diff forms]_i in new basis, GM is meant to preserve 4-dim subspace spanned by first 4 vectors.
		// We also calculate the relevant base change matrix Perm.
		Perm := ZeroMatrix(FieldOfFractions(R),d-1,d-1);
		for i := 1 to d-2 do
			if i lt ind then
				Perm[i,i] := 1;
			else
				Perm[ind,i] := -ResList[i+1]/ResList[ind];
				Perm[i+1,i] := 1;
			end if;
		end for;
		Perm[ind,d-1] := 1;
		M := Submatrix(Perm^(-1)*M*Perm,[1..d-2],[1..d-2]);
	end if;
	
	return M;
	
end function;

// This function returns the matrix of residues around a specified point. 

function ResidueMatrix(M: z:=0);

	R<x,t> := Parent(Numerator(M[1,1]));
	d := Degree(Parent(M));
	A := ZeroMatrix(Rationals(),d,d);
	for i := 1 to d do
		for j := 1 to d do
			M[i,j] := Evaluate(Numerator(M[i,j]),2,t+z)/Evaluate(Denominator(M[i,j]),2,t+z);
			res := Coefficient(Denominator(M[i,j]),2,1);
			if Coefficient(Denominator(M[i,j]),2,0) eq 0 and res ne 0 then
				A[i,j] := Coefficient(Numerator(M[i,j]),2,0)/res;
			end if;
		end for;
	end for;
	return A;
	
end function;
