/**********************************************************************************************
          POINT COUNTING ALGORITHM FOR CYCLIC COVERS OF THE PROJECTIVE LINE
            IMPLEMENTED BY CECILE GONCALVES - ECOLE POLYTECHNIQUE (LIX)
**********************************************************************************************/
/**********************************************************************************************
This program computes the Weil polynomial of a cyclic cover C : y^r = f(x) over a finite field, 
using Monsky-Washnitzer cohomology.
**********************************************************************************************/

forward pAdicIsZero;
forward pAdicDegree;
forward Iszero;
forward SimplifyCoeff;
forward Simplify;
forward Ser2Seq;
forward Normalize;
forward PrecompFrob;
forward Frob;
forward InvRthPower;
forward pAdicPolRingLiCombBaseCase;
forward pAdicPolRingLiComb;
forward FirstReduction;
forward SecondReduction;
forward ReductionCohomology;
forward MultiplyBlockMatrices;
forward FrobOfBlock;
forward NormOfBlockMatrix;
forward NormOfMatrix;
forward PermutationBlocks;
forward CharPolynomialFrobenius;
forward TestWeilPolynomialOK;
forward mult;


/**********************************************************************************************
 															CyclicCoverWeilPolynomial
**********************************************************************************************/
/* Returns the Weil polynomial of a cyclic cover C : y^r = f(x) of genus g over Fp^n 
The inputs are
r : The degree r of the cyclic cover
f : the defining polynomial of C, defined over Fp^n
Output : The numerator of the zeta function of C
*/
function CyclicCoverWeilPolynomial(r, f)
/* If f has a root in Fq, then we can reduce to the superelliptic case 
  Lroots := Roots(f);
  if (#Lroots ge 1) then
    alpha := Lroots[1][1];
    f1 := f div (Parent(f).1 - alpha);
	  h := Evaluate(f1, Parent(f).1 + alpha);
	  f := Parent(f)!(Parent(f).1^(Degree(f)-1)*Evaluate(h, 1/Parent(f).1));
	end if;*/

/* Step 0 : get the context : genus and basic rings */
  BaseField := CoefficientRing(f);
  p := Characteristic(BaseField);
  n := Degree(BaseField);
  P := DefiningPolynomial(BaseField);  
  d := Degree(f);
  delta := Gcd(r, d);
  g := ((r-1)*(d-1) - delta + 1) div 2;
  
 // Choose which basis we want to compute with : if we compute with B, then start = 2, otherwise start = 3
  if (p ge 2*r) then start := 3;
  elif (d gt r) then start := 3;
  else start := 2;
  end if;
  
  //print "start : ", start;
  
/* We will compute the Weil polynomial of C, then check if the result is ok. If so, we will return it, but if not, we will take larger bounds for precision */
  a := 0;
  check := false;
  
  while not check do
/* Step 1 : Precision bounds */
   // Weil bounds
   N0 := Ceiling(Log(p, 2*Binomial(2*g, g)) + (n*g/2));
 
   // Taking into account the loss of precision
   if (start eq 2) then
	   lbound := N0 + Ceiling(Log(p, ((d*p*(r-1) + r)/delta))) + Maximum(Ceiling(Log(p, (2*g - delta + 2)/delta)), Ceiling(Log(p, r))) + a;
	 else //start = 3 and we are computing with B'
     lbound := N0 + Ceiling(Log(p, 2*r-1)) + a;
   end if;
   
   N := lbound;
   if (start eq 2) then
     while (N - Floor(Log(p, p*(r*N - 1) - r)) lt lbound) do N +:= 1; end while;
     M := p*(N) - 1;
	 else //start = 3 and we are computing with B'
     while (N - Floor(Log(p, p*r*(N+1) - 3*r)) lt lbound) do N +:= 1; end while;
     M := p*(N + 1) - 2;
   end if;

   /*print "padic N0 :", N0;
   print "p adic N ", N;
   print "y adic", M;*/

// Construction of the rings we are computing with
  Zp := pAdicRing(p, N);
  Zpt<t> := PolynomialRing(Zp);
  Plift := Zpt!P;
  Zq<w> := ext<Zp|Plift>;
  Qq := FieldOfFractions(Zq);
  Zqx<x> := PolynomialRing(Zq);
  Zqxy<tau> := PowerSeriesRing(Zqx, M+1); // Here, tau = y^{-r}

// Lift the equation of C
	flift := Zqx!Polynomial([Zq!Eltseq(c): c in Coefficients(f)]);
	dflift := Derivative(flift); // derivative of flift

/* Step 2 : compute Frob(y)^-1 */
/* Precomputations */
// Lift the Frobenius to Zq
//Precompute Frob(w)
	wsigma := PrecompFrob(Plift, p, N, w);
//Compute fsigma = Frob(flift)
  L := Eltseq(flift);
  for i := 1 to (Degree(flift)+1) by 1 do 
    L[i] := Frob(L[i], p, wsigma);
    end for;
  fsigma := Polynomial(L);
  delete L;
  fsigma := Evaluate(fsigma, x^p);

// Compute R = 1 + (Frob(flift) - flift^p)*tau^p
  t3 := Cputime();
  R := Zqxy ! (1 + ((fsigma - flift^p)*(tau^p)));
  R := Normalize(R, flift, start);
//Newton iteration to get S = R^{-1/r}
  S := InvRthPower(R, r, N, flift, start);
  
  t4 := Cputime();
  //print "time for frob(y)^-1", t4-t3;

/* Step 3 : Action of Frobenius on cohomology */
/* Initialize the matrix of p-th power Frobenius to Zero Matrix : we will fill the good blocks later */
 	//MatrixOfFrobenius := ZeroMatrix(Zq, (r-1)*(d-1), (r-1)*(d-1));
 	MatrixOfFrobenius := ZeroMatrix(Qq, (r-1)*(d-1), (r-1)*(d-1));
/* Precomputations : compute PowersOfRoot such that PowersOfRoot[j] = Frob(y)^{-j}*/
  PowersOfRoot := [Zqxy | 0 : i in [1..r-1]];
  PowersOfRoot[1] := S;
  for j := 2 to r-1 do
  	PowersOfRoot[j] := Normalize(PowersOfRoot[j-1]*PowersOfRoot[1], flift, start);
  end for;
  delete R;
	delete S;
  
  // Precompute U and V such that Uflift + V dflift = 1, this is usefull for the Reduction Algorithm
  U, V := pAdicPolRingLiComb((Zqx ! 1), flift, Plift, p, N);

// for each vector in B or B' 	
  for j in [(r*(start - 2) + 1)..(r*(start - 2) + r-1)] do
    // For each j (for each subspace H^1(C)^{zeta_r^-j} of H^1(C) as in the article)
  	l := j*p mod r; // the image is in the subspace H^1(C)^{zeta_r^-l}
  	invfroby := Normalize(tau^(j*p div r)*(PowersOfRoot[1]^(r * (j div r))*PowersOfRoot[j mod r]), flift, start);
  	
	  for k in [0..d-2] do
	  // Compute the image of Frobenius (as a power series) of the vector x^k dx / y^{epsilon r + j}
  	  om := Ser2Seq(invfroby);
  	  for i in [1..#om] do 
  	    om[i] := om[i]*(p * x^(p*k + p - 1));
  	  end for;
  	  omega := Zqxy!om;
  	  omega := Normalize(omega, flift, start);
  	  om := Ser2Seq(omega);
  	  
  	 // Reduce the image into a linear combination of vectors of B (or B')
			T, NegativePowerOfp := ReductionCohomology(om, r, l, start, p, flift, U, V, N);

// Fill MatrixOfFrobenius
      for i := 0 to d-2 do
	      //MatrixOfFrobenius[(d-1)*(j mod r - 1) + (k + 1), (d-1)*(l-1) + i + 1] := Zq![(Integers()!co)/NegativePowerOfp : co in Eltseq(Coefficient(T, i))];
	      MatrixOfFrobenius[(d-1)*(j mod r - 1) + (k + 1), (d-1)*(l-1) + i + 1] := Qq![(Integers()!co)/NegativePowerOfp : co in Eltseq(Coefficient(T, i))];
	    end for;
  	end for;
  end for;
  
	  t5 := Cputime();
  	delete omega;
	  delete om;
	  delete invfroby;
	  delete PowersOfRoot;
	  //print "temps pour la réduction : ", t5-t4;

/* Step 4 : Compute the characteristic polynomial of the matrix */
// Compute the matrix of q-th power Frobenius, using the specific block struture of the matrix
	M, LblocksM, LposM := NormOfBlockMatrix(MatrixOfFrobenius, p, n, wsigma, d-1, [j*p mod r : j in [1..r-1]]);
	delete MatrixOfFrobenius;
// Compute The characteristic polynomial of M, using again the specific structure (block structure) of the matrix
	ZT<T> := PolynomialRing(Integers());
  CharPoly := ChangePrecision(CharPolynomialFrobenius(LblocksM, LposM, p^n, r), N0);

/* Step 5: Remove the extra Factor */
  if (delta gt 1) then
  // Compute the extra factor
	U := 1; 
		for i in [2..delta] do 
			if (delta mod i eq 0) then 
				a := 1; 
				while (p^(n*a) mod delta ne 1) do 
					a +:= 1;
				end while; 
				
				if (start eq 2) then
					U := U * (T^a - p^n)^(EulerPhi(i) div a); 
				else
					U := U * (T^a - 1)^(EulerPhi(i) div a); 
				end if;
			end if; 
		end for; 
		// Remove the extra Factor
		CharPoly := CharPoly div U;
		end if;

// Check that the coeff are ok
	CharPolynomial := ZT![(Integers()!coeff)mod p^N0: coeff in Eltseq(CharPoly)];
  ResultCoeff := [Integers() | 0 : i in [0..2*g]];
  ResultCoeff[1] := p^(n*g);
  ResultCoeff[2*g + 1] := 1;
   
   for i := 1 to g do
   		ci := Coefficient(CharPolynomial, 2*g - i) mod p^N0;
      if (ci lt ((p^(N0))/2)) then
         ResultCoeff[2*g - i+1] := ci;
      else
         ResultCoeff[2*g - i+1] := ci - p^N0;
      end if;
      ResultCoeff[i + 1] := p^(n*(g-i))*ResultCoeff[2*g - i+1];
   end for;

//Check that the Weil polynomial is OK. If it is, return the result ; otherwise take larger precision bounds
 check := TestWeilPolynomialOK(r, f, ReciprocalPolynomial(Polynomial(ResultCoeff)));
// if the test is not ok, we have to enlarge the precision bounds : 
a +:= 10;
end while;
	return ReciprocalPolynomial(Polynomial(ResultCoeff));
end function;

/**********************************************************************************************
 															Functions for p-adic manipulations
We know the precision bounds in theory, so we can make all the computations modulo p^N, where
N denote this precision and the computations will be done up to this precision, at the end. 
In Magma, as soon as we divide by p, we lose some precision digits, which make us lose some 
digits of precision in the computations. In order to avoid this phenomenon, we reajust the 
precision here.
**********************************************************************************************/
/* Returns whether a p-adic element over Z_q is 0 or not ( ie O(p^N))  */
pAdicIsZero := function(coeff)
	L := Eltseq(coeff);
	res := (Rationals() ! L[1]) eq (Rationals() ! 0);
	i := 2;
	
	while (res and i lt #L + 1) do
	  res := (Rationals() ! L[i]) eq (Rationals() ! 0);
		i := i + 1;
	end while;
	
	return res;
end function;

/* Returns the degree of a polynomial whose coefficients are in Zq (the O(p^stuff) are considered as 0*/
pAdicDegree := function(pol)
  d := Degree(pol);
  if(d ne -1) then 
    while (d gt -1 and pAdicIsZero(Coefficient(pol, d))) do
  	  d := d-1;
	  end while;
  end if;
  return d;
end function;

/* Returns whether a polynomial over Zq is zero or not */
Iszero := function(pol)
  d := pAdicDegree(pol);
  return (d eq -1);
end function;

/* Simplify a p-adic element (erase all the O(p^stuff)) and hence avoids unfair loss of precision.
It returns the element coeff of Zq after simplification. */
SimplifyCoeff := function(coeff)
  L := Eltseq(coeff);
  for i in [1..#L] do 
	  L[i] := Universe(L)!Rationals()!L[i];
	end for;
  return Parent(coeff)!L;
end function;

/* Simplify a polynomial over Zq (uses SimplifyCoeff).
It returns the univariate polynomial pol defined over Zq after simplification */
Simplify := function(pol)
  L := Eltseq(pol);
  L := L[1..pAdicDegree(pol)+1];

  if (Iszero(pol)) then return Parent(pol)!0;
  else
    for i in [1..#L] do 
    	L[i] := SimplifyCoeff(L[i]);
    end for;
  return Polynomial(L);
  end if;
end function;

/**********************************************************************************************
 															Function for manipulate power series
**********************************************************************************************/
/* Convert a power series into a list.
Returns a list L such that L[i] is the (i+1)-th coeff of the power series ser in input. */
Ser2Seq := function(ser)
  d := Degree(ser);
  L := [Parent(Coefficient(ser, d))|];
  for i in [0..d] do
  	L := L cat [(Coefficient(ser, i))];
  end for;
  return L;
end function;

/**********************************************************************************************
         									Reduction Functions and tools for reduction
These functions are useful for reducing an element of the first group of cohomology represented
by a power series into a linear combination of basis vectors
**********************************************************************************************/
/* Uses the equation of the curve to normalize an element of Adagger, that is, to reduce the terms of 
degree in x greater than Degree(f).
It returns a power series which is equivalent (in the cohomology) to the input ser 
Here, the input flift is the defining polynomial of the curve and start indicates which set of 
differentials of the cohomology we are computing with (see article, Section 7) */
Normalize := function(ser, flift, start)
  d := Degree(flift);
  L := Ser2Seq(ser);

  // Use the equation of the curve at each coefficient to reduce the degree in x when it is greater than the degree of flift
  for i:= #L to start by -1 do
    Bi := L[i];
    if (Degree(Bi) gt d-1) then
      Qi, Ri := Quotrem(Bi, flift);
      L[i] := Ri;
      L[i-1] := (L[i-1] + Qi);
    end if;
  end for;
  return Parent(ser)!L;
end function;

/* These two following functions compute polynomials over Zq such that hlift = Uflift + Vdflift, which is used for Step 1  in Algorithm 2 (reduction algorithm) 
when we need to compute A and B such that Aflift+Bflift' = R.
This is done by using the extended euclidean algorithm in Fq[x] and by computing the coefficients of U and V for each power of p:
We can see U(x) (respectively V(x)) as U(x) = \sum_{i = 0}^{n-1}(U_i(x) p^i)) where the Ui(x) are polynomials over Fq.
*/
/* Computes U0 and V0 such that U0*f + V0 f' = h mod p (over Fq). U0 and V0 aue to precision 1.
The inputs are:
hlift, flift: Polynomials such that Uflift + V dflift = hlift
Plift : The defining polynomial of the unramified extension Qq of Qp
p : the characteristic of Fq
*/
pAdicPolRingLiCombBaseCase := function(hlift, flift, Plift, p)
   Fp := FiniteField(p);
   Fpt<t> := PolynomialRing(Fp);
   P := Fpt ! Plift;
   Fq<w> := ext<Fp | P>;
   Fqx<x> := PolynomialRing(Fq);
   f := Fqx ! [Fq ! [Integers()!co : co in Coefficients(coeff)] : coeff in Coefficients(flift)];
   df := Fqx ! [Fq ! [Integers()!co : co in Coefficients(coeff)] : coeff in Eltseq(Derivative(flift))];
   h := Fqx ! [Fq ! [Integers()!co : co in Coefficients(coeff)] : coeff in Eltseq(hlift)];
   g, l, k := XGCD(f, df);
   quot, rem := Quotrem(k*h, f);
   U_1 := Parent(flift) ! [CoefficientRing(flift) ! Eltseq(coeff) : coeff in Eltseq(l*h + quot*df)];
   V_1 := Parent(flift) ! [CoefficientRing(flift) ! Eltseq(coeff) : coeff in Eltseq(rem)];
   return U_1, V_1;
end function;      

/* Computes U and V such that U*flift + V0 dflift = hlift over Zq
The inputs are:
hlift, flift: Polynomials such that Uflift + V dflift = hlift
Plift : The defining polynomial of the unramified extension Qq of Qp
p : the characteristic of Fq
N: the p-adic precision
*/
pAdicPolRingLiComb := function(hlift, flift, Plift, p, N)
   U, V := pAdicPolRingLiCombBaseCase(hlift, flift, Plift, p);
   
   for i := 2 to N do
      res := (hlift - U*flift - V*Derivative(flift)) div p^(i-1);
      U_0, V_0 := pAdicPolRingLiCombBaseCase(res, flift, Plift, p);
      U := U + p^(i-1)*U_0;
      V := V + p^(i-1)*V_0;
   end for;
   
   return Parent(hlift) ! U, Parent(hlift) ! V;
end function;      

/* Step 1 of Algorithm 2 (ReductionCohomology Algorithm) : reduces the degree in tau
Given an element of the cohomology om = \sum{i le mu}{A_i(x) y^{-ri + \epsilon r + l} dx, it returns an element S(x) dx/y^{\epsilon r + l} equivalent to om 
Note that \epsilon = 0 if we are computing with B and \epsilon = 1 if we compute with B'
The inputs are:
om : the differential form we want to reduce
r: the exponent such that C: y^r = f(x)
l: the exponent such that \sum{i le mu}{A_i(x) y^{-ri + \epsilon r + l} dx
start : start = 2 if we are computing with B, and 3 if we use B'
p : the chartacteristic of Fp
flift : the defining polynomial of C (over Z_q)
U, V : are such that Uflift + V flift' = 1
*/
FirstReduction := function(om, r, l, start, p, flift, U, V)
		for k := #om to start by -1 do
			Rk := om[k];
			// Computes Ak and Bk such that Rk = Ak flift + Bk flift'
     	if (Degree(Rk) lt 2*Degree(flift) - 1) then 
 	    	Ak := (Rk * U) mod Derivative(flift);
   	  	Bk := (Rk * V) mod flift;
     	else
     		Ak := (Rk * U);
	    	Bk := (Rk * V);
	    end if;
	    // Replaces the term Rk tau^k dx/{y^{\epsilon*r + l}} with (Ak + (r/(r(k-1)+l)*Bk') tau^{k-1} dx/{y^{\epsilon*r + l}}
	    if ((Ak ne 0) or (Bk ne 0)) then 
	    	if (Valuation((r*(k-2)+l), p) eq 0) then 
	    		om[k-1] := om[k-1] + Simplify((Ak + (Universe(om)!(r/(r*(k-2)+l)))*Derivative(Bk)));
	   	  else 
	    		om[k-1] := om[k-1] + Simplify(Ak + (Universe(om)!Polynomial([CoefficientRing(flift)![((r*(Rationals()!coeff))/(r*(k-2)+l)) : coeff in Eltseq(co)] : co in Coefficients(Derivative(Bk))])));
	    	end if;
			end if;			    
	    om[k] := Universe(om)!0;
 	   end for;
 	   
 	   return om[start-1];
end function;

/* Step 2 of Algorithm 2 (ReductionCohomology Algorithm) : reduces the degree in x
Given an element of the cohomology of the form om = S(x) dx/y^{\epsilon r + l}, returns an element T(x) dx/y^{\epsilon r + l} equivalent to om, with Degree(T) less than d-1.
Note that \epsilon = 0 if we are computing with B and \epsilon = 1 if we compute with B'
The inputs are:
S : the polynomial such that the differential form we want to reduce is S(x) dx/y^{\epsilon r + l}
r: the exponent such that C: y^r = f(x)
l: the exponent such that \sum{i le mu}{A_i(x) y^{-ri + \epsilon r + l} dx
start : start = 2 if we are computing with B, and 3 if we use B'
flift : the defining polynomial of C (over Z_q)
p : the chartacteristic of Fp
N : The p-adic precision
*/
SecondReduction := function(S, r, l, flift, p, N)
		T := S;
    m := Degree(T);
    NegativePowerOfp := 1;
    d := Degree(flift);
	  g := ((r-1)*(d-1) - GCD(r, d) + 1) div 2;
    x := Parent(S).1;
    
    while (m gt d-2) do 
      // Computes Ttilde as in the article
      if (m eq d-1) then 
      	Tt := Derivative(flift);
      else
	      Tt := (r*(m - (d-1))*x^(m - (d-1) - 1)*flift + (r-l)*x^(m - (d-1))*Derivative(flift));
	    end if;
	    
	    // Puts the good coefficient in Ttilde in order its leading coefficient is the same as the leading coefficient of T and reduces the degree in x
	    if (Valuation((Coefficient(Tt, m))) eq 0) then 
	    	Tt := Parent(T)!((Coefficient(T, m))/(Coefficient(Tt, m)))*Tt; //Puts the good leading coefficient
	      T := Simplify(T - Tt); // Reduces the degree in x
	    else //Puts the good leading coefficient but with a storage of the denominators 
	    	LC := Coefficient(Tt, m);
	    	PowOfp := Gcd(Eltseq(LC) cat [p^N]);
	    	LC := LC div PowOfp;
	    	Tt := Parent(S)!(1/LC)*Tt;
	    	LCT := Coefficient(T, m);
	    	PowOfpT := GCD(GCD([Integers()!co : co in Eltseq(LCT)]), p^N);
	    	
	    	if (PowOfpT ge PowOfp) then 
	    		Multiply := CoefficientRing(T)!(LCT/PowOfp);
	    	else
	    		PowOfDiv := PowOfp div PowOfpT;
	    		NegativePowerOfp *:= PowOfDiv;
	    		T := Parent(T)!(PowOfDiv)*T;
	    		Multiply := CoefficientRing(T)![(Integers()!co) div PowOfp : co in Eltseq(Coefficient(T, m))];
	    	end if;
	    	
	    	T := Simplify(T - (Parent(T)!Multiply)*Tt); // Reduces the degree in x   	
	    end if;
	    	    
      m := Degree(T);
    end while;
    
    return T, NegativePowerOfp;
end function;

/* Algorithm 2 (ReductionCohomology Algorithm) : reduces a differential form into a linear combination of vectors of B or B'
Given an element of the cohomology of the form om = S(x) dx/y^{\epsilon r + l}, returns an element T(x) dx/y^{\epsilon r + l} equivalent to om, with Degree(T) less than d-1.
Note that \epsilon = 0 if we are computing with B and \epsilon = 1 if we compute with B'
The inputs are:
S : the polynomial such that the differential form we want to reduce is S(x) dx/y^{\epsilon r + l}
r: the exponent such that C: y^r = f(x)
l: the exponent such that \sum{i le mu}{A_i(x) y^{-ri + \epsilon r + l} dx
start : start = 2 if we are computing with B, and 3 if we use B'
flift : the defining polynomial of C (over Z_q)
p : the chartacteristic of Fp
N : The p-adic precision
*/
ReductionCohomology := function(om, r, l, start, p, flift, U, V, N)
// Step 1
	F := FirstReduction(om, r, l, start, p, flift, U, V);
// Step 2
  F, NegativePowerOfp := SecondReduction(F, r, l, flift, p, N);
  return F, NegativePowerOfp;
end function;

/**********************************************************************************************
         															Lifting functions
These functions are usefull to lift the Frobenius to Zq and cohomology.
**********************************************************************************************/
/* Lift of Frobenius map to p-adics elements:
An element of Zq is seen as a degree less than n polynomial over Zp: z = sum{z_k w^k}.
We will first precompute a lift of Frob(w) and then compute Frob(z) by Horner's method (see Section 5 of the article)
 */
/* This function returns the lift of w by a Newton iteration, because w is a root of P, the defining polynomial
of the unramified extention Qq of Qp:
x_{n+1} := x_n - P(x_n)/P'(x_n)
It returns Frob(w), the lift of the Frobenius map applied on w.
The inputs are:
P : the defining polynomial of the unramified extention Qq of Qp
p: the characteristic of Fp
prec: the p-adic precision required
w: elements of Zq are seen as elements of Zp[w] modulo P
*/
PrecompFrob := function(P, p, prec, w)
  // Compute the Frobenius on w modulo p
	res := Parent(w)![(Integers()!co) mod p : co in Eltseq(w^p)];
	
	// Apply a Newton iteration to find a root of P modulo the required precision
	for i in [1..Ceiling(Log(2, prec))] do
		res := res - ((Evaluate(P,res))/(Evaluate(Derivative(P),res)));
	end for;
	
	return res;
end function;

/* This function computes the lift of Frobenius applied to an element of Zq by Horner's method
It returns Frob(z), the lift of the Frobenius map applied on z, an element of Zq.
The inputs are:
elt : the element of Zq we want to apply the Frobenius : elt = \sum{elt_k w^k}
p: the characteristic of Fp
wsigma: The Frobenius applied to w
*/
Frob := function(elt, p, wsigma)
	L := Eltseq(elt);
  res := L[#L];
  
  for i := #L - 1 to 1 by -1 do
  	res := res*wsigma + L[i];
  end for;
  
  return res;
end function;

/* Lift of Frobenius map to cohomology:
An element of the cohomology is seen as a degree less than M polynomial over Zq[x] (it it an overconvergent power series whose coefficients are
univariate polynomials over Zq, up to precision N ; as the series is overconvergent, at some point all the coefficients are O(p^N), that is 0 in Zq).

We will first precompute a lift of Frob(y^-1) and then compute Frob(omega) where omega is a "basis vector".
 */
/* This function returns the inverse of the r-th squareroot of elt (an element of Adagger) given in input by a Newton iteration:
It computes a root of P(X) := elt*X^r - 1 by a Newton iteration : 
x_{n+1} := x_n - P(x_n)/P'(x_n) = x_n - (elt*x_n^r - 1)/(r*elt*x_n^{r-1})
As we want a root of P, we can replace elt*x_n^{r-1} by x_n^{-1}
So the iteration is: x_{n+1} := x_n - x_n(elt*x_n^r - 1)/r
The inputs are:
elt : the element we want the inverse of a r-th squareroot
r: the power of the squareroot
prec: the p-adic precision required
flift: the defining polynomial of the curve (for normalisation) of cohomological elements
start: indicate which set of differential is used for the computations : if start=2, then we use the canonical set of differentials B, else, if start = 3, then we use B' (see Section 7 of the Article)
*/
InvRthPower := function(elt, r, prec, flift, start)
  S := Parent(elt) ! 1;
  for i:=1 to Ceiling(Log(2, prec)) by 1 do 
  	T := Normalize(S^r, flift, start); //x_n^r
  	T := Normalize(elt*T, flift, start); // elt*x_n^r
  	T := Normalize((r+1)-T, flift, start); // (r+1) - elt*x_n^r
  	T := (S*T)/r; // x_n*[(r+1) - elt*x_n^r]/r = x_n + x_n/r - x_n*elt*X_n^r/r = x_n - x_n(elt*x_n^r - 1)/r
    S := Normalize(T, flift, start);
  end for;
  return S;
end function;

/**********************************************************************************************
						Compute the characteristic polynomial using the block structure of the matrix
**********************************************************************************************/
/* Computes the "Norm" of the Matrix Mat in input, that is, M*Msigma*...*Msigma^n-1 which is
the matrix of the q-th power Frobenius when Mat is the matrix of the p-th power Frobenius.
This function works for a block squared matrix Mat which has exactly one nonzero squared block on 
each row and each column. 
The inputs are:
Mat : Matrix we want to compute the norm (here this is the matrix of the p-th power Frobenius)
p : characteristic of Fp
n : degree of the extension Fq of Fp
wsigma : the frobenius of w, the generator of Zq as an unramified extension of Zp
sizeblock : size of the squared blocks
Lposblock: List of the positions of the non-zero block on each column : the non zero block on the j-th column
is the sizeblock*sizeblock squared matrix on the block line Lposblock[j]
*/
NormOfBlockMatrix := function(Mat, p, n, wsigma, sizeblock, Lposblock)
	LRes := [Submatrix(Mat, (Lposblock[j] - 1 )*sizeblock + 1 , (j-1)*sizeblock + 1, sizeblock, sizeblock) : j in [1..#Lposblock]]; //Excludeing the blocks to apply Frobenius on it
	LMatFrob := [FrobOfBlock(LRes[i], p, n, wsigma) : i in [1..#LRes]]; // Frobenius on blocks
	LposRes := Lposblock;

  // starting multiplying
  for i in [1..n - 1] do 
    LRes, LposRes := MultiplyBlockMatrices(LMatFrob, Lposblock, LRes, LposRes);
  end for;
  
  // Get the result
  MatrixResult := ZeroMatrix(CoefficientRing(Mat), Nrows(Mat), Ncols(Mat));
  for j in [1..#LposRes] do
    InsertBlock(~MatrixResult, LRes[j], (LposRes[j] - 1)*sizeblock + 1, (j - 1)*sizeblock + 1);
  end for;
 
  return MatrixResult, LRes, LposRes;
end function;

/* Given a matrix Block, returns the matrix whose coefficients are the Frobenius applied on the coefficients of FrobOfblock */
FrobOfBlock := function(Block, p, n, wsigma)
	L := Eltseq(Block);
	dim := Nrows(Block);
	
	  	for j in [1..#L] do 
      	L[j] := Frob(L[j], p, wsigma);
    	end for;
  return Matrix(Universe(L), dim, dim, L);
end function;

/* Computes the product of the two block matrices A and B, which are both block matrices with exactly one nonzero block on each
row and each column
The inputs are:
LblocksA : The list of the non zero blocks of A : LblocksA[j] is the block on the j-th block column
LposA: The list of the position of non zero blocks of A : LposA[j] is the block line index of the nonzero block at block column j
LblocksB : The list of the non zero blocks of B
LposB: The list of the position of non zero blocks of B : LposB[j] is the block line index of the nonzero block at block column j
The result is a matrix C = AB returned under the form LblocksC, LposC
*/
MultiplyBlockMatrices := function(LblocksA, LposA, LblocksB, LposB)
  LblocksC := [Universe(LblocksA)|];
  LposC := [];
  
  for j in [1..#LposA] do
    // We have C_{i, j} = A_{i, k}*B_{k, j} where k is the index such that B_{k, j} is the non-zero block
    // Find the block index i of the block line where the nonzero block on block colums j lies
    i := LposA[LposB[j]];
    LposC := LposC cat [i];
    LblocksC := LblocksC cat [LblocksA[LposB[j]]*LblocksB[j]]; // multiply the right blocks
  end for;
  
  return LblocksC, LposC;
end function;

/* Computes the characteristic polynomial of the q-th power Frobenius:
The matrix of the q-th power Frobenius is a block matrix with exactly one nonzero block on each row and column.
We can reorganize its rows block and columns blocks such that we have a block diagonal matrix with small blocks in each diagonal
block corresponding to the permutation j*q mod r.
So we proceed into 3 steps:
1. Reorganize the matrix regarding to the permutation j*q mod r : we will get a list of blocks
2. Compute the characteristic polynomial of each big block
3. The characteristic polynomial of the Matrix of Frobenius is the product of the latter characteristic polynomials.

The inputs are
LblocksM: list of nonzero blocks of M such that LblocksM[j] is the nonzero block on column j
LposM: LposM[j] is the index of the block line of the nonzero block on the j-th column
q: the order of Frobenius map (M is the matrix of q-th power Frobenius
r: the order of cyclic cover

It returns the characteristic polynomial of M given under the form LblocksM, LposM
*/
CharPolynomialFrobenius := function(LblocksM, LposM, q, r)
  sizeblock := Nrows(LblocksM[1]);
/* Step 1 : Reorganize the lines and columns of M (it doesn't change the characteristic polynomial) with respect to the permutation of block colums j -> jq mod r */
// Compute the permutation
	Lindex, Llengths := PermutationBlocks(q, r);
// Compute the corresponding big diagonal blocks
  Lblocks := []; //list of big diagonal blocks
  currentindex := 1; //start of the cycle we are dealing with
  for j in [1..#Llengths] do //for each big diagonal block (composed of little blocks)
    Block := ZeroMatrix(CoefficientRing(LblocksM[1]), Llengths[j]*sizeblock, Llengths[j]*sizeblock); //initialize the big block matrix
    for i in [1..Llengths[j] - 1] do //fill the matrix with respect to cycles
      InsertBlock(~Block, LblocksM[Lindex[currentindex]], (i)*sizeblock + 1 , (i-1)*sizeblock + 1);
      currentindex +:= 1;
    end for;
    InsertBlock(~Block, LblocksM[Lindex[currentindex]], 1 , (Llengths[j]-1)*sizeblock + 1);
    currentindex +:= 1;
    Lblocks := Lblocks cat [Block];
  end for;
  
/* Step 2: compute the characteristic polynomials of the diagonal blocks */
  LCharpoly := [CharacteristicPolynomial(B) : B in Lblocks];

/* Step 3: The characteristic polynomial of the Matrix Of Frobenius is the product of the previous characteristic polynomials */
  CharPoly := LCharpoly[1];
  for i in [2..#LCharpoly] do
    CharPoly *:= LCharpoly[i];
  end for;
  
  return CharPoly;
end function;

/* Returns Lindex, Llengths which are respectivly the list on index and lenghts of the cycles of the 
permutation of {1, ... , r*1} : j -> j*q mod r

q: the order of Frobenius map (M is the matrix of q-th power Frobenius
r: the order of cyclic cover
*/
PermutationBlocks := function(q, r)
  E := {i : i in [1..r-1]};
  Lindex := [];
  Llenghts := [];
  
  while (#E gt 0) do
	  e := Minimum(E); //Initializes the beginning of a new cycle
	  Exclude(~E, e);
	  Lindex := Lindex cat [e]; // Add the index to the list of indexes
	  c := 1; //initializes the size of the new cycle
	  
	  i := e*q mod r;
	  Exclude(~E, i);
		
		while (i ne e) do
			Lindex := Lindex cat [i]; // add i to the list of indexes
			c +:= 1; //	increase the length of the cycle
			i := i*q mod r; //continue the cycle
	 		Exclude(~E, i);
	 	end while;
	 	Llenghts := Llenghts cat [c]; //Add the new size of a cycle (because we just finished a cycle)
	 end while;
	 
   return Lindex, Llenghts;
 end function;

/* Matrix of q-th power frobenius, normal case (ie without using the block structure) which is not used here */
NormOfMatrix := function(Mat, p, n, wsigma)
	Mt := Mat;
	L := Eltseq(Mat);
	dim := Nrows(Mat);
	
	  	for j in [1..#L] do 
      L[j] := Frob(L[j], p, wsigma);
    end for;
    M := Matrix(Universe(L), dim, dim, L);
	
  for i in [1..n - 1] do 
    Mt := M * Mt;
  end for;
  
  return Mt;
end function;

/**********************************************************************************************
						Check that the given polynomial is the good one, that is, corresponds to the
 														Weil polynomial of the curve
**********************************************************************************************/
/* 
Check that the polynomial WeilPoly given in input is the Weil Polynomial of the curve C : y^r = f(x)
defined over F_p^n
Inputs : 
r : the index r in the defining equation of C
f : the defining polynomial of C over F_p^n
WeilPoly : the polynomial we want to check

Output : returns true if it is ok, false otherwise

*/
TestWeilPolynomialOK := function(r, f, WeilPoly)
	// Get the hypothetic cardinality of the Jacobian of C: y^y = f
	N := Evaluate(WeilPoly, 1);
	
	// Define the spaces and the curve
	Fq := CoefficientRing(Parent(f));
	A<X, Y> := AffineSpace(Fq, 2);
	C := Curve(A, Y^r - Evaluate(f, X));

  // Construct a degree 0 divisor (not principal)
  isprincipal := true;
  while (isprincipal) do
  	// take some random points
	  p1 := RandomPlace(C, 1);
		p2 := RandomPlace(C, 1);
		p3 := RandomPlace(C, 1);
		//Construct the degree 0 divisor
		D := (p1 + p2 - 2*p3);
		isprincipal := IsPrincipal(Reduction(D));
	end while;
	
	//Compute N*D and check that the result is principal	
	K := mult(N, D);
	
	return IsPrincipal(K);
end function;

/* Given an integer N and a divisor D, computes N*D as a reduced divisor
Inputs : 
N and D
Output : The reduced divisor N*D
*/
mult := function(N, D)
	if N eq 1 then return Reduction(D);
	else
		if (N mod 2 eq 0) then return $$(N div 2, Reduction(2*D));
		else return Reduction(D + $$(N div 2,Reduction(2*D)));
		end if;
	end if;
end function;


