{kedlaya(Qpolydach,p,fdach,debug)=

   /* 
      This program computes the rational points of a hyperelliptic curve y^2 = Q(x)
      following a algorithm of Kedlaya[2001].


      Input:
      Qpolydach: The seperable polynomial Q(x) from F_q[x] with y^2 = Q(x).
      p: The characteristic of F_q.
      n: The power of p, ie. q = p^n.
      fdach: The polynomial f(x) from F_p[x] with F_q ~ F_p[x] / <f(x)>.
      debug: The debug-modus, 0 is no debugging up to ???

      Output:
      The number of rational points of the hyperelliptic curve.


      Dieses Programm berechnet die rationalen Punkte auf einer
      hyperelliptischen Kurve y^2 = Q(x) nach dem Algorithmus 
      von Kedlaya[2001]. 

      Eingabe:
      Qpolydach: (Polynom) Das seperable Polynom Q(x) aus F_q[x] mit y^2 = Q(x).
      p: (Primzahl) Die Charakteristik von F_q.
      n: (positiver Integer) Der Exponent von p, dh.q =  p^n.
      fdach: (Polynom) Das Polynom f(x) aus F_p[x] mit F_q ~ F_p[x] / <f(x)>.
      debug: (Boolean??) Der Debugger-Modus, von 0 (kein debugging) bis ??.

      Ausgabe:
      Die Anzahl der rationalen Punkte der hyperelliptischen Kurve
      in F_q.
   */
   local(Qpoly,mod,modulou,Qpolydachv,vorlaeufigpraezisionN,ggT,SUM,vek,prez,frobprez);
  
   y;
   Y;
   w;  
   v;
   u;
   /* Die benoetigte Ordnung der Variablen fuer die GP-Funktion Mod */

    /* 
      Ich veraendere die Variablen, weil x in PARI eine fuer mich ungeeignete Sonderstellung hat...
      Also fdach(z) -> fdach(u) und Qpolydach(z,x) -> Qpolydach(u,v).  
    */
    
    print(" ");
    print("----------------------------------------------------------------------------------------------------------------");
    print("in kedlaya.gp  ");
    print("----------------------------------------------------------------------------------------------------------------");
    
    modulou = subst( fdach , variable( fdach ) , u );
    
    Qpolydachv = subst( subst( lift(Qpolydach) , variable( lift(Qpolydach) ) , v) , variable(fdach) , u )* Mod( 1 , modulou );

    ggT = bezout(Qpolydachv, deriv(Qpolydachv,v));

    /*  - Hier untersuchen wir, ob Qpolydach und fdach geeignet sind. */
    
    if(
       poldegree(lift(ggT[3]),u) == 0,
       /* 
         Falls der Grad vom ggT Null ist, so ist der ggT eine Konstante, also invertierbar und somit 1.
       */
       
       
       if(
          issquarefree( modulou ),
	  if( 
             length(factormod( modulou , p )) != 2,
	     /*  Falls mehrere Faktoren auftauchen  */
             error( "Das erzeugende Polynom f(x)= "fdach" ist nicht irreduzibel! " ),
            ),
	     /* Falls doppelte Faktoren auftauchen  */
	    error("Das erzeugende Polynom f(x)= "fdach" ist nicht irreduzibel! ");  /*  endif   */
            );  /*  endif  */
       print("Polynom in Ordung, Erzeuger auch"),
       error("Polynom Q(x) hat mehrfache Nullstellen");
      ); /*  endif  */  

   vorlaeufigpraezisionN = ceil( ( poldegree( Qpolydachv  , v ) ) * log(2) / log(p) + poldegree( modulou , u )*( poldegree( Qpolydachv  , v ) - 1 ) / 4); 

   prez = vorlaeufigpraezisionN + valuation( poldegree( Qpolydachv  , v ) , p ) 
   + ceil( log( poldegree( Qpolydachv , v ) - (2 / p)) / log(p) ) + ceil(log(poldegree( Qpolydachv , v ) - 2)/log(p));
   
   frobprez = 0;  

   /* 
      Da die Praezision > 0 und floor(0) = 0.  
   */
   
   while(
        frobprez -floor( log( 2*frobprez + 1 ) / log(p) ) < prez,
        frobprez++
         );  /* endwhile */
frobprez = frobprez + debug;


/*  Lifte die Daten nach Zp bzw. Zq  */

   mod = lift(modulou)+O(p^prez);

   Qpoly = Mod(1+O(p^prez),mod)*(lift(lift(Qpolydachv)));


vektor=frob(ggT[2],Qpoly,mod,p,prez,frobprez);

matrize(vektor,Qpoly,mod,p,prez)   

}




{addhelp(kedlaya,"kedlaya(Qpol,p,modulo,{b}) counts the points of the hyperelliptic curve y^2=Q(x) "
                 "over F_q ~ F_p[z]/<modulo(z)> using the algorithm of Kedlaya[2001]. Therefore the prime "
                 "has to be odd. It returns a vector with the number of points in the first entry "
                 "and in the second the Zeta-function of the curve."); 
}

{yinverse(Qpoly,modulo,p,prez,frobprez)=


/*  

   This part computes the frobeniusaction on 1/y. Here we use y instead of 1/y.

   Input:
   
   Qpoly = The lifted polynomial in Zq[v] with y^2 = Qpoly.
   modulo = Tthe lifted polynomial in Zp[u] with Qq ~ Qp[u]/(modulo(u)).
   p = The characteristic of the residue field of Qq.
   prez = The necessary precision for the computation in Qq.
   frobprez = The necessary number of summands in the frobeniusexpansion.

   Output:
   
   inv = The Polynomial with the frobeniusaction on y.

   Dieser Teil berechnet die Frobeniusaktion auf 1/y. Wir benutzen y anstatt 1/y.

   Eingabe:
   
   Qpoly = Das geliftete Polynom in Zq[v] mit y^2 = Qpoly.
   modulo = Das geliftete Polynom in Zp[u] mit Qq ~ Qp[u]/(modulo(u)).
   p = Die Charakteristik des Restklassenkoerpers Fq.
   prez = Die benoetigte Praezision fuer die Rechnungen in Qq.
   frobprez = Die Anzahl der benoetigten Summanden in der Frobeniusaktion.

   

Hier invertiere ich y^2 von y^2=Q(x).

Hier benutze ich zum ersten mal die Variable y! Y := y^(-1)
*/

local(inv,zwischen,zwischen2,E,zwischen3,testzwischen,bi,pow);




print("----------------------------------------------------------------------------------------");
print("In yinverse");
print("----------------------------------------------------------------------------------------");

gettime();

E = smallfrob(Qpoly,u+O(p^prez),modulo,p,prez,1);

inv = Mod(1+O(p^prez),modulo);
 
zwischen = substpol(E,Qpoly,x^2);

zwischen = Mod(1+O(p^prez),modulo) + subst( zwischen,x,Y)*y^(2*p);

zwischen = sorting(zwischen,Qpoly,modulo,p,prez,frobprez);

bi = length(binary(2*p*frobprez+1));

forstep(
        b = bi-1,
        0,
        -1,
         print("In yinverse:: In forstep-Schleife "b" bis 0");
                           
         pow = ceil((2*p*frobprez+1)/(1<<b));
 
         inv = truncate((3/2)*inv-(1/2)*zwischen*inv^3 + O(y^pow));

         /*  Die Newton Iteration fuer die inverse Quadratwurzel  */
     
         inv = subst(substpol(inv,Qpoly,x^2),x,Y);

         /*  Das Reduzieren der Koeffizienten bezueglich des Grades  */
         
         inv = sorting(inv,Qpoly,modulo,p,prez,frobprez);

         /*  Sortiere die entstandenen Y ein  */
 
         print("Dauer "gettime()" und die Prez "padicprec(inv,p)" und der Grad in y "poldegree(inv,y));
    );  /*  endforstep  */

inv = lift(lift(inv))*Mod(1+O(p^prez),modulo);

return(inv*y^p);

}





{frob(mu,Qpoly,modulo,p,prez,frobprez)=
	
/*
  This part of the program computes the action of the p-frobenius on the basis.

  Input:

  mu = The polynomial in Fq[v] with 1 = lambda Qpoly + mu Qpoly' mod p.
  Qpoly = The lifted polynomial in Zq[v] with y^2 = Qpoly.
  modulo = Tthe lifted polynomial in Zp[u] with Qq ~ Qp[u]/(modulo(u)).
  p = The characteristic of the residue field of Qq.
  prez = The necessary precision for the computation in Qq.
  frobprez = The necessary number of summands in the frobeniusexpansion.
  
  Output:

  frobdx = the vector with the reduced basis (reduced by the part "reduction"), where
           the first entry is the frobeniusaction on x^0(dx/y).

  Dieser Teil des Programms berechnet die Aktion des p-Frobenius auf der Basis.

  Eingabe:
  mu = Ein polynom aus Fq[v] mit 1 = lambda Qpoly + mu Qpoly' mod p.
  Qpoly = Das geliftete Polynom in Zq[v] mit y^2 = Qpoly.
  modulo = Das geliftete Polynom in Zp[u] mit Qq ~ Qp[u]/(modulo(u)).
  p = Die Charakteristik des Restklassenkoerpers Fq.
  prez = Die benoetigte Praezision fuer die Rechnungen in Qq.
  frobprez = Die Anzahl der benoetigten Summanden in der Frobeniusaktion.

  Output:
  
  frobdx = Ein Vektor mit der reduzierten Basis (reduziert von reduction), also der erste Koeffizient ist die Anwendung
           auf das Basiselement x^0(dx/y) und so weiter.

*/
local(ggt,frobdx,zwischen);

ggt = liftbezout(Qpoly,mod,mu,p,prez);

yinv = yinverse(Qpoly,mod,p,prez,frobprez);

print("Berechne jetzt die Frobeniusaktion");

frobdx = vector(poldegree(Qpoly,v)-1);

for(
    n = 0,
    poldegree(Qpoly,v)-2,
    
    print("in der "n"'ten for Schleife von "poldegree(Qpoly,v)-2" in der Berechnung der Frobeniusaktion");
    
    zwischen = subst(substpol(p*v^(p*n+p-1)*yinv,Qpoly,x^2),x,Y);  
    
    /*  Anwendung des Frobenius auf das Basiselement x^n(dx/y) und die Reduzierung der Grade der Koeffizienten */
    
    zwischen = Mod(1+O(p^prez),mod)*lift(Mod(1,p^prez)*lift(lift(zwischen)));

    /* Schmeisse ueberfluessige Praezision weg  */
    
    zwischen = sorting(zwischen,Qpoly,mod,p,prez,frobprez);
 
    /* Sortiere die negativen Potenzen von y, also Y ein  */
  
    zwischen = reduction(zwischen,ggt,Qpoly,mod,p,prez);
    
    /* Reduziere die Aktion zurueck auf die Basis  */

    frobdx[n+1] = zwischen;
    
    );  /*  endfor  */


kill(yinv);

print("++++++++++++++++++++++++++++++++++++++++++++++");
print("verlasse frobreduction");
print("++++++++++++++++++++++++++++++++++++++++++++++");

return(frobdx);

}

{reduction(ser,ggt,Qpoly,mod,p,prez)=

/*
  This part reduces the frobeniusaction back to the basis.

  Input:
  
  ser = The Polynomial in Qq[v,y,Y], which is the action of the frobenius on a basis element.
  ggt = A vektor with the liftet Bezoutcoefficients.
  Qpoly = The lifted polynomial in Zq[v] with y^2 = Qpoly.
  mod = The lifted polynomial in Zp[u] with Zq ~ Zp[u]/(mod(u)).
  p = The characterisic of the residue field Fq of Zq.
  prez = The nessecary precision for the computation in Zq.
 
  Output:

  SER = The Polynomial in Z[v], which is the reduction to the basis.

*/
local(SER,Koeff,poldiv,dQpoly,gradmesser);

print("------------------------------------------------------------");
print("in reduction");
print("------------------------------------------------------------");

dQpoly = deriv(Qpoly,v);


SER = polcoeff(ser,poldegree(lift(lift(ser))),y);


forstep(
           j = poldegree(lift(lift(ser)),y)-2,
	   1,
	   -2,
   
	   poldiv1 = divrem(simplify(SER*ggt[2]),Qpoly);
	
           if(
              type(poldiv1[1])=="t_RFRAC",
              print("divrem hat wieder das Ergebnis nicht berechnet");
              print("Versuche es erneut!");
              print("  ");
              poldiv1 = divrem(simplify(SER*ggt[2]),simplify(Qpoly));
              if(
                 type(poldiv1[1])=="t_RFRAC",
                 error("divrem will nicht!");
                );
              );
           
           /* Da divrem manchmal anstatt mit Rest zu teilen, einfach die rationale Funktion SER*ggt[2]/Qpoly ausgibt, fange ich hier den Fehler ab */
	   
            Koeff = SER*ggt[1]+ poldiv1[1]*dQpoly + (2/(j))*deriv(poldiv1[2],v);
    
	    SER = Koeff + polcoeff(ser,j,y);
	     
	     
	     );  /*  endforstep  */



/* Es koennen noch negative Potenzen in Form von Y in SER sein  */

if(
   poldegree(lift(lift(ser)),Y) > 0,
   
   forstep(
              h = poldegree(lift(lift(ser)),Y),
	      1,
	      -2,
	      
	      if(lift(lift(polcoeff(ser,h,Y)))==0,

              /*  Falls der Koeffizient Null ist, tue nichts  */
	         ,
             
	      SER = SER + (polcoeff(ser,h,Y)*Qpoly^((h-1)/2 + 1));
	     
              /*  Ansosnsten hebe den Koeffizienten an  */
 
	      );  /*  endif  */

	);  /*  endforstep  */

   
);  /*  endif  */



SER = lift(lift(SER))*Mod(1+O(p^(padicprec(SER,p))),modulo);

forstep(
           k = poldegree(lift(lift(SER)),v),
           poldegree(Qpoly,v),
	   -1,
	   
	  
	  Koeff =2*deriv(v^(k+1-poldegree(Qpoly,v)),v)*Qpoly + v^(k+1 -poldegree(Qpoly,v))*dQpoly; 
	  
	  SER = SER - ((polcoeff(SER,k,v)/(pollead(Koeff,v)))*(2*deriv(v^(k+1-poldegree(Qpoly,v)),v)*Qpoly + v^(k+1 -poldegree(Qpoly,v))*dQpoly));
	 
          /*  Hier reduziere ich die Grade in v  */
 
	);  /*  endforstep  */

if(
   poldegree(lift(lift(SER)),v) >= poldegree(Qpoly,v),

   gradmesser = poldegree(dQpoly,v),

   gradmesser = poldegree(lift(lift(SER)),v);

  );


/* Weil poldegree manchmal die Grd nicht richtig zurueck gibt, setze ich ihn hier notgedrungen */

if(
   gradmesser==poldegree(dQpoly,v),
   
    
   SER = SER - (polcoeff(SER,poldegree(dQpoly),v)/poldegree(Qpoly,v))*dQpoly;

/*  Hier benutze ich dy = Qpoly'(dx/2y)  */ 
  );
  
  /* Nun ist SER auf die Basis reduziert und wird somit zurueck gegeben  */
  
print("Verlasse reduktion  ");

return(SER);
   
   
}

{sorting(pol,Qpoly,modulo,p,prez,frobprez)=

/*  
  This part sorts the Y, which is 1/y, in a Polynomial in y.

  Input:

  pol = A polynomial in Qq[v,y,Y].
  Qpoly = The lifted polynomial in Zq[v] with y^2 = Qpoly.
  modulo = Tthe lifted polynomial in Zp[u] with Qq ~ Qp[u]/(modulo(u)).
  p = The characteristic of the residue field of Qq.
  prez = The necessary precision for the computation in Qq.
  frobprez = The necessary number of summands in the frobeniusexpansion.
  
  Output:

  zwischen = The sorted polynomial.

  Dieser Teil des Programms sortiert die Y in das Polynom pol ein, dabei ist Y eigentlich 1/y.

  Eingabe:
  
  pol = Ein Polynom in Qq[v,y,Y].
  Qpoly = Das geliftete Polynom in Zq[v] mit y^2 = Qpoly.
  modulo = Das geliftete Polynom in Zp[u] mit Qq ~ Qp[u]/(modulo(u)).
  p = Die Charakteristik des Restklassenkoerpers Fq.
  prez = Die benoetigte Praezision fuer die Rechnungen in Qq.
  frobprez = Die Anzahl der benoetigten Summanden in der Frobeniusaktion.

  Output:
  
  zwischen = Das sortierte Polynom.
  
*/

local(zwischen);




print("----------------------------------------------------------------------------------------");
print("In sorting");
print("----------------------------------------------------------------------------------------");

gettime();


zwischen = 0;
for(z = 0,
        poldegree(lift(lift(pol)),y),
       
	if(
	    poldegree(polcoeff(lift(lift(pol)),z,y),Y) > 0,
	  
	    if(
	        poldegree(polcoeff(lift(lift(pol)),z,y),Y) <= z,

		for(
	            q = 0,
		    poldegree(polcoeff(lift(lift(pol)),z,y),Y),
		    zwischen = zwischen + polcoeff(polcoeff(pol,z,y),q,Y)*y^(z-q);
		    ),  /* endfor  */
		
               
		for(
		    q = 0,
		    z,
		    zwischen = zwischen + polcoeff(polcoeff(pol,z,y),q,Y)*y^(z-q);
		    );  /* endfor  */
		for(
		     q = z+1,
		     poldegree(polcoeff(pol,z,y),Y),
		     zwischen = zwischen + polcoeff(polcoeff(pol,z,y),q,Y)*Y^(q-z);
		     );  /*  endfor  */
		     
		
	       ),  /*  endif  */
	       
	       if(
	          poldegree(polcoeff(pol,z,y),Y) == 0,
	         
	          zwischen = zwischen + polcoeff(pol,z,y)*y^z;
		  );  /*  endif  */
	       
	       
	       
	       );  /*  endif  */
	       
   );  /*  endfor   */

return(zwischen);

}

{smallfrob(Qpoly,zero,modulo,p,prez,bool)=

/*
   This part computes the frobenius on Qq or the E= sigma_p(Qpoly)-Q^p.

   Input:

   Qpoly = The lifted polynomial in Zq[v] with y^2 = Qpoly.
   zero = A generator of the extension Qq/Qp.
   modulo = Tthe lifted polynomial in Zp[u] with Qq ~ Qp[u]/(modulo(u)).
   p = The characteristic of the residue field of Qq.
   prez = The necessary precision for the computation in Qq.
   bool = A boolean variable. 
  
  Output:

  if bool = 0, then it gives zeromodp back, which is the the generator of Qq/Qp with zeromodp = zero^p mod p.
  if bool = 1, then it gives E back. 


  Dieser Teil des Programms berechnet den Frobenius auf Qq oder das E = sigma_p(Qpoly)-Q^p.

  Eingabe:
  
  Qpoly = Das geliftete Polynom in Zq[v] mit y^2 = Qpoly.
  zero = Ein Erzeuger der Erweiterung Qq/Qp. 
  modulo = Das geliftete Polynom in Zp[u] mit Qq ~ Qp[u]/(modulo(u)).
  p = Die Charakteristik des Restklassenkoerpers Fq.
  prez = Die benoetigte Praezision fuer die Rechnungen in Qq.

  Output:
  
  Falls bool = 0, so gibt das Programm zeromodp zurueck, welches ein Erzeuger von Qq/Qp mit zeromodp = zero^p mod p.
   

*/


local(summe,zeromodp,E,sigmaQpoly,coeff,Qhochp,dmodulo,inverse,bi,testost);

dmodulo = deriv( modulo , u );
bi = binary(prez);
liftmodulo = lift(modulo);
zeromodp = Mod(lift(lift(Mod(1,p)*Mod(lift(lift(zero))^p, modulo))),liftmodulo);

/*  Die Initialisierung von der Nullstelle   */

inverse = Mod(lift(lift(Mod(1,p)*subst(lift(dmodulo),u,zeromodp)^(-1))),liftmodulo);

/* Die Initialisierung von modulo(zero)'^(-1)   */

forstep(
        i= length(bi)-1,
	0,
	-1,
	pow = ceil((prez)/(1<<i));
	
		
        /* die Newton Iteration vom Frobenius */
	zeromodp = Mod(1,p^pow)*(zeromodp - subst(modulo,u,zeromodp)*inverse);
		
		
	/*  Die Newton Iteration fuer die Inverse */
	inverse = Mod(1,p^pow)*(2*inverse - subst(dmodulo,u,zeromodp)*inverse^2);
		    
	zeromodp = Mod(lift(lift(zeromodp)),modulo);
	inverse = Mod(lift(lift(inverse)),modulo);
                   
       );  /*  endforstep  */

zeromodp = (lift(lift(zeromodp)))*Mod(1+O(p^prez),modulo);

if(
   bool == 0,
   return( zeromodp ),
   

   sigmaQpoly = subst(subst(lift(Qpoly),u,zeromodp),v,v^p);


   E = sigmaQpoly - Qpoly^p;
   return(E);
  );  /*  endif  */    
}
   
{liftbezout(Qpoly,mod,lambda,p,prez)=

/*

  This part computes the bezout-coefficients of gcd(Qpoly,deriv(Qpoly,v)) via
  a Newton iteration from the gcd mod p.
  
  Input:
  
  Qpoly = The polynomial in Zq[v] with y^2 = Qpoly.
  lambda = the polynomial in Fq[v]  which is a  the bezout-coefficient of the
           Reduction of the gcd. Hence Qpolymodp*mu + deriv(Qpolymop,v)*lambda = gcdmodp.
  mod = The polynomial in Zp[u] with Zq ~ Zp[u]/<mod>
  p = The characteristic of the residue field Fq of Zq
  prez = The necessary precision for the computation in Zq.

  Output:
  
  [Mu,Lambda] = the vector with the bezout-coefficients of the gcd(Qpoly,deriv(Qpoly,v)), hence
                Mu*Qpoly + Lambda*deriv(Qpoly,v)= gcd.

  
  Dieser Teil berechnet die Bezoutkoeffizienten des ggT(Qpoly,Qpoly') mittels
  einer Newton-Iteration von dem ggT mod p.

  Eingabe:

  Qpoly = Das geliftete Polynom in Zq[v] mit y^2 = Qpoly.
  lambda = Das Polynom in Fq[v], das der Bezoutkoeffizient vor Qpolymodp ist.
  mod = Das geliftete Polynom in Zp[u] mit Qq ~ Qp[u]/(mod(u)).
  p = Die Charakteristik des Restklassenkoerpers Fq.
  prez = Die benoetigte Praezision fuer die Rechnungen in Qq.

  Output:
   
  [Mu,Lambda] = Der Vektor mit den gelifteten Bezoutkoeffizienten des ggTs.


*/
local(modmodp,inv,grad,liftlambda,liftQpoly,dliftQpoly,liftprez);



liftlambda = lift(lift(lambda));  /* liftlambda ist ein Polynom in 2 Variablen aus Z[v]  */



liftQpoly = lift(lift(Qpoly));   
dliftQpoly = deriv(liftQpoly,v);
dQpoly = deriv(Qpoly,v);
liftmod = lift(mod);
   
liftprez = length(binary(prez))-1;
inv = [0,Mod(1,mod)*(liftlambda+O(p^prez))];

forstep(
        i= liftprez,
	0,
	-1,
	print("In liftbezout:: In Newton Iteration Nr. "i);
	powprez = ceil(liftprez/(1<<i));
        inv = divrem( 2*inv[2] - dQpoly*inv[2]^2,Qpoly);     
        /* die eigentliche Newton Iteration  */
 
       );  /*  endforstep  */
   
   grad = poldegree(inv[2],v);
  
   Lambda = (inv[2]+O(p^prez))*Mod(1,mod);

   qr = divrem(Mod(1+O(p^prez),mod)-Lambda*deriv(Qpoly,v),Qpoly);

   if(
      (poldegree(qr[2] , v )) < 0,
   
      Mu = qr[1],
      error("Falsche Division in liftbezout "qr[2]);
     );  /*  endif  */

   return([Mu,Lambda]);

}


{matrize(vek,Qpoly,mod,p,prez,pot)=

/*
  This part computes the characteristc polynomial fo the frobenius.

  Input:
  
  vek = A vector with the reduced frobeniusaction on the basis.


*/

local(mat,Xpol,frobZq,zwischenmat,var1,var2,pow);

print("----------------------------------------------------------------------------------------------------");
print("Berechne Matrix und die Zeta-Funktion");
print("----------------------------------------------------------------------------------------------------");

mat= matrix(poldegree(Qpoly,v)-1,poldegree(Qpoly,v)-1,X,Y,polcoeff(vek[Y],X-1,v));


frobZq = smallfrob(Qpoly,u+O(p^prez),mod,p,prez,0);

zwischenmat = mat * (subst(lift(mat),u,frobZq)); 

for(
    j = 2,
    poldegree(mod,u)-1,                        
    frobZq = smallfrob(Qpoly,frobZq,mod,p,prez,0);
    zwischenmat = zwischenmat*(subst(lift(mat),u,frobZq)); /*  die Aktion des Frobenius  !! */ 
   );  /*  endfor  */

Xpol = charpoly(zwischenmat,x);

print("Das charpoly ungeliftet "Xpol);
print("  ");
print("Das charpoly!! geliftet "lift(Mod(1,p^prez)*lift(lift(Xpol))));


Ypol=0;

for(
    i=1,
    (poldegree(Qpoly,v)-1)/2-1,

    print("Inner for-Schleife: "i);
    
    pow = ceil(log(2*binomial(poldegree(Qpoly,v)-1,i))/log(p) + (i/2)*poldegree(mod,u));
    
    print("Wat is denn die Potenz pow: "pow);

    var1 = lift(Mod(1,p^pow)*lift(lift(polcoeff(Xpol,poldegree(Qpoly,v)-1-i,x))));

if(
       type(var1)=="t_POL",
       var1 = polcoeff(var1,0,v);
       if(
          type(var1)=="t_POL",
          var1 = polcoeff(var1,0,u);
          if(
             type(var1)=="t_POL",
             var1 = polcoeff(var1,0,y);
             if(
                type(var1)=="t_POL",
                var1 = polcoeff(var1,0,Y);
               );
             );

         );  /*  endif  */

     );  /*  endif  */
    
   
    var2 = ceil((p^pow)/2);   

    if(
       var2 < var1, 
       print("Musste das Vorzeichen anpassen");
       var1 = var1-p^pow;
      );  /*  endif  */


    Ypol = Ypol + var1*y^(poldegree(Qpoly,v)-1-i)+ 
              var1*p^(poldegree(mod,u)*((poldegree(Qpoly,v)-1)/2-i))*y^i;
 
  );  /*  endfor  */


pow = ceil(log(2*binomial(poldegree(Qpoly,v)-1,(poldegree(Qpoly,v)-1)/2))/log(p) + ((poldegree(Qpoly,v)-1)/4)*poldegree(mod,u));    




var1 =lift(Mod(1,p^pow)*lift(lift(polcoeff(Xpol,((poldegree(Qpoly,v)-1)/2),x))));



if(
       type(var1)=="t_POL",
       var1 = polcoeff(var1,0,v);
       if(
          type(var1)=="t_POL",
          var1 = polcoeff(var1,0,u);
          if(
             type(var1)=="t_POL",
             var1 = polcoeff(var1,0,y);
             if(
                type(var1)=="t_POL",
                var1 = polcoeff(var1,0,Y);
               );
             );

         );  /*  endif  */

     );  /*  endif  */


var2 = ceil((p^pow)/2);

print("var1 :"var1);


if(
   var1 > var2,
   var1 = var1 - p^pow;
  
   if(debug==1,
       print("Musste das Vorzeichen anpassen");
       ),
);  /*  endif  */
print("var1 :"var1);

Ypol = Ypol + y^(poldegree(Qpoly,v)-1) + var1*y^((poldegree(Qpoly,v)-1)/2) + p^(poldegree(mod,u)*((poldegree(Qpoly)-1)/2));

if(
   abs(polcoeff(Ypol,1,y)) > p^(poldegree(mod,u)),
   
   var1 = p^(poldegree(mod,u)) + polcoeff(Ypol,poldegree(Qpoly)-2,y) + 1,
   
   var = p^(poldegree(mod,u)) + polcoeff(Ypol,1) + 1;
  );  

/* Da polcoeff in meinen Beispielen den Koeffizienten hier (und nur hier...)  nicht immer (eigentlich nie) richtig ausgibt, ist diese seltsame Abfrage notwendig!  */


return([polrecip(Ypol),var1]);


}