/* unitcomp.m -- joint computation with Romyar */

intrinsic Eta(p::RngIntElt, z::.) -> .
{}
   r := 2;
   while Numerator(BernoulliNumber(r)) mod p ne 0 and r lt p+100 do
      r := r+2;
   end while;
   if r ge p then
      r := p-3;  // for testing.
   end if;
   print "Making eta (r=",r,")";
   eta := &*[(&+[z^j : j in [0..i-1]])^(i^(p-1-r) mod p) : i in [1..p-1]];
   //eta := eta0 * (1-z)^(&+[i^(p-1-r) : i in [1..p-1]] mod p);
   return eta;
end intrinsic;

intrinsic RelativeExtension(p::RngIntElt) -> FldNum
{}
   // construct Qzetap
   Qzetap<z> := CyclotomicField(p);
   // construct eta
   eta := Eta(p,z);
   S<y> := PolynomialRing(Qzetap);
   print "Making relative extension";
   time L<a> := quo<S|y^p-eta>;   // L = Qzetap(eta^(1/p)).
   return L;
end intrinsic;

intrinsic LocalEmbeddings(p::RngIntElt, L::RngUPolRes) -> .
{}
   R<x> := PolynomialRing(RationalField());
   f := Evaluate(R!CyclotomicPolynomial(p),x+1);   // Eisenstein polynomial
   Qzetap_p<t> := TotallyRamifiedExtension(pAdicField(p), f) ;   // t = zeta_p - 1.
   g := Modulus(L);
   S<X> := PolynomialRing(Qzetap_p);
   phi := hom<CyclotomicField(p) -> Qzetap_p | t+1>;      // zeta_p |---> t+1 = (zeta_p - 1) + 1 = zeta_p.
   phi_g := &+[phi(Coefficient(g,n))*X^n : n in [0..Degree(g)]];
   print "Now factor defining polynomial over Qp(zeta_p)...";
   // time F := Factorization(phi_g);
   psi := hom<CyclotomicField(p) -> GF(p) | 1>;      // zeta_p |---> 1 (mod p)
   P<Y> := PolynomialRing(GF(p));
   phi_g_modp := &+[psi(Coefficient(g,n))*Y^n : n in [0..Degree(g)]];
   print "phi_g_modp = ", phi_g_modp;
   return phi_g, phi_g_modp;
   F := Factorization(phi_g_modp);

   /* In the p=37 example, the roots all have valuation 16, according to MAGMA's ValuationsOfRoots command. */

end intrinsic;


intrinsic FindRoot(p::RngIntElt) -> .
{}
   f := CyclotomicPolynomial(p);
   O := EquationOrder(f);
   z := O.2;
   I := ideal<O | 1-z>;      // maximal ideal of O.
   Q, phi := quo<O|I^2>;
   zbar := phi(z);
   print "making eta";
   time eta := Eta(p,zbar);
   print "computing elements of O/I^2";
   time Qset :=[a+b*zbar : a in [0..36], b in [0..36]];
   assert #Set(Qset) eq #Q;
   return eta, Set([x : x in Qset | x^p eq eta]);
end intrinsic;



///////////////////////////////////////////

intrinsic FindDegree37Field() -> SeqEnum
{}
   r := 32;
   p := 37;
   del_pow := 2;

   function res(n)
      return n mod p;
   end function;
   function respow(r,s) 
      return Integers()!((IntegerRing(p)!r)^s);
   end function;

   K<z> := CyclotomicField(p);
   eta := &*[(1-z^i)^respow(i,p-r-1) : i in [1..p-1]];
   R<X> := PolynomialRing(K);
   print "Forming field L.";
   time L<w> := quo<R | X^p - eta>;
   print "L = ", L;

   function del(x)
      assert x in L;
      del_of_w := (w^respow(2,r)) * 
        &*[(1-z^k)^((res(respow(k,p-r-1)*respow(del_pow,r))-respow(k,p-r-1)*respow(del_pow,r)) div p) : k in [1..p-1]];
      phi := hom<L->L | del_of_w>;
      y := phi(x);
      elt := Eltseq(y);
      del_of_elt := [Conjugate(e,z^del_pow) : e in elt];
      return L!del_of_elt;
   end function;

   return L, del;

   function random_K()
      print "Constructing random element of K...";
      return K![Random(-1,1) : i in [1..Degree(K)]];
   end function;

   function random_L()
      print "Constructing random element of L...";
      return L![random_K() : i in [1..p]];
   end function;

   alpha := random_L();

   print "Computing trace of random element (alpha =", alpha,")";
   t := del(alpha);
   tr_alpha := alpha + t;
   for i in [1..p-3] do
      printf "i = %o (of %o).\n",i, p-3;
      time t := del(t);               // t = del^(i+1)(alpha).  
      tr_alpha := tr_alpha + t;
   end for;
   print "Finished computing trace of random element.";   

   print "Computing minimal polynomial of trace(alpha), as an element of K[X]";
   time f := MinimalPolynomial(tr_alpha);
   return f;

   // compute charpoly of trace.  This will be a poly of degree p*(p-1) that should be a perfect (p-1)st power.



   /* The following fails with 
      i = 30 (of 36).
      i = 31 (of 36).
      malloc block past in more_data scan
      Magma: Internal error
      Time to this point: 320.67
    */

   assert false;

   print "Computing another random element gamma.";
   gamma := random_L();
   print "gamma = ", gamma;
   print "Computing images of gamma under powers of tr(alpha).";
   t := tr_alpha*gamma;
   ims := [gamma, t];   
   for i in [1..p-1] do
      printf "i = %o (of %o).\n",i, p-1;
      time t := tr_alpha*t;   // t = tr(alpha)^(i+1)*gamma.
      Append(~ims, t);
   end for;
   print "Finished computing images.";
   
   ////////////////////////////////////////////////////////////
   // find a linear relation between the elements of ims.
   ////////////////////////////////////////////////////////////
   print "Finding a linear relation between the images";
   V := VectorSpace(RationalField(),p*(p-1));
   function vectorize(x)   // view element of L as vector over Q.
      return V![Coefficient(Coefficient(x,i),j) : i in [0..p-1], j in [0..p-2]];
   end function;

   print "(finding matrix)";
   ims_vec := [vectorize(x) : x in ims];   // list of 37 vectors in a VS of dimension 37*36 over Q.
   A := Transpose(MatrixAlgebra(p,p*(p-1))!&cat[Eltseq(v) : v in ims_vec]);
   print "(computing kernel)";   
   time ker := Kernel(A);
   print "Kernel = ", ker;  // this should have dimension 1.
   return Basis(ker);

end intrinsic;




intrinsic FindDegree37Field0() -> SeqEnum
{}
   r := 32;
   p := 37;
   del_pow := 2;

   function res(n)
      return n mod p;
   end function;
   function respow(r,s) 
      return Integers()!((IntegerRing(p)!r)^s);
   end function;

   K<z> := CyclotomicField(p);
   eta := &*[(1-z^i)^respow(i,p-r-1) : i in [1..p-1]];
   R<X> := PolynomialRing(K);
   print "Forming extension L.";
   time L<w> := quo<R | X^p - eta>;
   print "L = ", L;
   function MyConjugate(x, zpow)
      psi := hom<K -> K | zpow>;
      return psi(x);
   end function;

   del_of_w := (w^respow(del_pow,r)) *
              &*[(1-z^k)^((res(respow(k,p-r-1)*respow(del_pow,r))
               -respow(k,p-r-1)*respow(del_pow,r)) div p) : k in [1..p-1]];
   function del(x)
      assert x in L;
      ans := L!0;
      del_of_w_pow := 1;
      for n in [0..p-1] do 
         an := Coefficient(x,n);
         if an ne 0 then
            ans := ans + MyConjugate(an,z^del_pow)*del_of_w_pow;
         end if;
         del_of_w_pow := del_of_w_pow*del_of_w;
      end for;
      return ans;
   end function;

   function random_L()
      print "Constructing 'random' element of L...";
      return L![(-1)^i : i in [1..p]];
   end function;

   alpha := w;

   function trace(x)
      printf "%o steps: ", p-3;
      t := del(x);
      tr_x := x + t;
      for i in [1..p-3] do
         printf "%o, ",i, p-3;
         t := del(t);               // t = del^(i+1)(x).  
         tr_x := tr_x + t;
      end for;
      return tr_x;
   end function;

   print "Computing trace of random element (alpha =", alpha,")";
   time tr_alpha := trace(alpha);
   print "\nFinished computing trace of random element.";   

   print "Computing minimal polynomial of trace(alpha), as an element of K[X]";
   time f := MinimalPolynomial(tr_alpha);

   return f, L, K, del, tr_alpha;
end intrinsic;

///////////////////////////////////////////

intrinsic FindDegree37FieldMod(ell::RngIntElt) -> SeqEnum
{}
   r := 32;
   p := 37;
   del_pow := 2;


   function res(n)
      return n mod p;
   end function;
   function respow(r,s) 
      return Integers()!((IntegerRing(p)!r)^s);
   end function;

   R_ell := PolynomialRing(GF(ell));
   cyclo := R_ell!CyclotomicPolynomial(p);
   require Order(Integers(p)!ell) eq 36  : "Argument 1 must be inert in Q(zeta_p).";

   K<z> := quo<R_ell | cyclo>;
   eta := &*[(1-z^i)^respow(i,p-r-1) : i in [1..p-1]];
   R<X> := PolynomialRing(K);
   print "Forming extension L.";
   time L<w> := quo<R | X^p - eta>;
   print "L = ", L;
   function MyConjugate(x, zpow)
      psi := hom<K -> K | zpow>;
      return psi(x);
   end function;

   del_of_w := (w^respow(del_pow,r)) *
              &*[(1-z^k)^((res(respow(k,p-r-1)*respow(del_pow,r))
               -respow(k,p-r-1)*respow(del_pow,r)) div p) : k in [1..p-1]];
   function del(x)
      assert x in L;
      ans := L!0;
      del_of_w_pow := 1;
      for n in [0..p-1] do 
         an := Coefficient(x,n);
         if an ne 0 then
            ans := ans + MyConjugate(an,z^del_pow)*del_of_w_pow;
         end if;
         del_of_w_pow := del_of_w_pow*del_of_w;
      end for;
      return ans;
   end function;


   function random_L()
      print "Constructing 'random' element of L...";
      return L![(-1)^i : i in [1..p]];
   end function;

   alpha := w;

   function trace(x)
      printf "%o steps: ", p-3;
      t := del(x);
      tr_x := x + t;
      for i in [1..p-3] do
         printf "%o, ",i, p-3;
         t := del(t);               // t = del^(i+1)(x).  
         tr_x := tr_x + t;
      end for;
      return tr_x;
   end function;

   print "Computing trace of random element (alpha =", alpha,")";
   time tr_alpha := trace(alpha);
   print "\nFinished computing trace of random element.";   

   print "Computing minimal polynomial of trace(alpha), as an element of K[X]";
   time f := MinimalPolynomial(tr_alpha);
   print "f = ", f;
   S<T> := PolynomialRing(GF(ell));
   g := S!0;
   for i in [1..38] do
       assert Degree(Coefficient(f,i-1)) le 0;
       g := g + T^(i-1)*Coefficient(Coefficient(f,i-1),0);
   end for;
   print "g = ", g;
   return g, L, K, del, tr_alpha;

end intrinsic; 


forward CRT_vec, LiftToZ;

intrinsic Find37FieldCRT(stop::RngIntElt) -> RngUPolElt
{}
   // Warning -- remove ``Proof = false'' to make this rigorous.
   M := [* *];
   polcyclo := CyclotomicPolynomial(37);
   
   tm := Cputime();
   t := 0;
   ell := 13; 
   while true and ell le stop do
      while not IsIrreducible(PolynomialRing(GF(ell))!polcyclo) do 
         ell := NextPrime(ell); 
      end while;      
      f := FindDegree37FieldMod(ell);
      if Degree(f) eq 37 then
         Append(~M,Eltseq(f));  
         t +:= 1;
         if t mod 5 eq 0 then
            v := CRT_vec(M); 
            vlift := LiftToZ(v) ;
            if vlift eq LiftToZ(M[1]) then
               break;
            end if;
            M := [* v *];
            print "CRT step -- ", Cputime(tm);
            print "So far:  ", vlift;
            tm := Cputime();
         end if;
      end if;
      ell := NextPrime(ell);
   end while;
   print "Found it!!!";
   return PolynomialRing(Integers())!LiftToZ(v);
end intrinsic;

function LiftToZ(T)
   if Type(T) eq AlgMatElt then
      return MatrixAlgebra(Rationals(),Degree(Parent(T)))!LiftToZ(Eltseq(T));
   end if;
   assert Type(T) eq SeqEnum;

   d := #T;
   N := Characteristic(Parent(T[1]));

   T0 := [];
   for i in [1..d] do
      if (Integers()!T[i]) ge N/2 then     // or something like this.
         T0[i] := (Integers()!T[i])-N;
      else
         T0[i] := Integers()!T[i];
      end if;
   end for;

   return T0;
end function;

function CRT_vec(list)
   // input is a list of vectors of the same degree modulo various coprime numbers
   assert Type(list) eq List;
   assert #list gt 0;

   d := #list[1];
   assert d gt 0;
   M := [Characteristic(Parent(v[1])) : v in list];
   return [Integers(LCM(M))| CRT([Integers()!v[i] : v in list], M) : i in [1..d] ];
end function;

