Q := Rationals();
Z := IntegerRing();
Qt<T> := PolynomialRing(Q);
Zt<U> := PolynomialRing(Z);

intrinsic respow(p:: RngIntElt, a :: RngIntElt, b :: RngIntElt) -> RngIntElt 
{return the residue mod p of a^b}
  return Z!((GF(p)!a)^b);
end intrinsic;

intrinsic KummerExt(p :: RngIntElt, r :: RngIntElt) -> .
{returns data on the Kummer extension of Q(zeta_p) assoc. to r}

  cyclo := Qt!CyclotomicPolynomial(p);
  K<z> := quo<Qt | cyclo>;
  eta := &*[(1-z^i)^respow(p,i,p-r-1) : i in [1..p-1]];
  Kt<X> := PolynomialRing(K);
  L<w> := quo<Kt | X^p - eta>;
  prim := PrimitiveRoot(p);
  
  function del(x)
     del_of_w := w^respow(p,prim,r) *  
     &*[(1-z^k)^(((respow(p,k,p-r-1)*respow(p,prim,r) mod p)
     -respow(p,k,p-r-1)*respow(p,prim,r)) div p) : k in [1..p-1]];
     elt := Eltseq(x);
     del_of_elt := &+[hom<K->K|z^prim>(elt[k]) * del_of_w^(k-1) : k in [1..#elt]];
  
     return L!del_of_elt;
  end function;
  
  return w,del;
end intrinsic;


intrinsic TraceElt(p :: RngIntElt, r :: RngIntElt, w :: RngUPolResElt, del) 
-> RngUPolResElt
{returns an element of a degree p extension F of L (hopefully a generator)}
  
  d := del(w);
  tr := w + d;
  for i in [1..p-3] do
      d := del(d);               // d = del^(i+1)(x).  
      tr +:= d;
  end for;

  return tr;
end intrinsic;

intrinsic MinPoly(p :: RngIntElt, r :: RngIntElt, tr :: RngUPolResElt) -> RngUPolElt
{returns a minimal polynomial over Q of the element tr}
   prim := PrimitiveRoot(p);
 
   function tracemod(l)
	Rl := PolynomialRing(GF(l));
	cyclol := Rl!CyclotomicPolynomial(p);
	Kl<zl> := quo<Rl | cyclol>;
	etal := &*[(1-zl^i)^respow(p,i,p-r-1) : i in [1..p-1]];
	Sl<Xl> := PolynomialRing(Kl);
	Ll<wl> := quo<Sl | Xl^p - etal>;
	
	function dell(x)
  		del_of_wl := wl^respow(p,prim,r) * 
  		   &*[(1-zl^k)^(((respow(p,k,p-r-1)*respow(p,prim,r) mod p)
     		-respow(p,k,p-r-1)*respow(p,prim,r)) div p) : k in [1..p-1]];
  		elt := Eltseq(x);
  		del_of_elt := &+[hom<Kl->Kl|zl^prim>(elt[k]) * del_of_wl^(k-1) : k in [1..#elt]];
  		return Ll!del_of_elt;
	end function;

	function tracel(x)
     	   t := dell(x);
     	   tr_x := x + t;
     	   for i in [1..p-3] do
        	t := dell(t);               // t = del^(i+1)(x).  
        	tr_x +:= t;
     	   end for;
     	   return tr_x;
	end function;
	
	return tracel(wl);
   end function;

   function minpoly()
	l := NextPrime(5000000);
	done := false;
	N := 1;
	g := 0;
	h := 0;

	function getcoeff(f,i)
		a := Eltseq(Coefficient(f,i));
		if #a eq 1 then return Z![a[1]];
		else return 0;
		end if;
	end function;
		
	function smabs(a,M)
		if (a le M/2) then return a;
		else return a-M;
		end if;
	end function;	
		
	while (not done) do
		if IsPrimitive(IntegerRing(p)!l) then 
			trl := tracemod(l);
			if (#Eltseq(trl) gt 1) then
				f := MinimalPolynomial(trl);
				g := h;
				M := N*l;
				if (g eq 0) then
					h := &+[ smabs(getcoeff(f,i),M)*U^i : i in [0..p]];
				else
					h := &+[ smabs(CRT([getcoeff(f,i), Z!Coefficient(g,i)],[l,N]),M)*U^i 
					: i in [0..p]];
				end if;

				// h;
				// print "\n";
				if (g eq h) then 
					// print "equal";
					done := (Evaluate(g,tr) eq 0); 
				end if;
				N := M;
			end if;
		end if;
		l := NextPrime(l);
	end while;
	
	return g;
   end function;

   f := minpoly();
   return f;
end intrinsic;

intrinsic ImprovePoly(p :: RngIntElt, f :: RngUPolElt) -> .
{replaces f(x) by g(x) = f(x/m) for m maximal such that f has integer coefficients}
   Re := RealField();
   tot := Floor(Minimum([Root(AbsoluteValue(Coefficient(f,i)),p-i) : i in [0..p-1] 
	  | not (Coefficient(f,i) eq 0) ]));
   m := 1;
   a := 1;
   while (a le tot) do
      a +:= m;
      if {Coefficient(f,i) mod a^(p-i) : i in [0..p-1]} eq {0} then 
          m := a;
      end if;
   end while;   
   
   g := &+[(Coefficient(f,i) div m^(p-i))*U^i : i in [0..p]];
   return g,m;
end intrinsic; 

intrinsic MaxOrder(p :: RngIntElt, g :: RngUPolElt) -> RngOrd
{computes the maximal order in F}
   
   E := EquationOrder(g);
   Ep := pMaximalOrder(E, p);
   i := Abs(Z!Discriminant(Ep));
   i := i div p^Valuation(i,p);
   _, i := IsSquare(i);
   
   T := TrialDivision(i,10000000);
   m := Maximum([T[j][2] : j in [1..#T]]);
   if m ge 2 then
   		for n in [j : j in [1..#T] | T[j][2] ge 2] do
                     Ep := pMaximalOrder(Ep, T[n][1]);                     
                end for;
   		i := Abs(Z!Discriminant(Ep));
                i := i div p^Valuation(i,p);
   		_, i := IsSquare(i);
   end if;

   I := pRadical(E, i);          // this computes the "p" maximal
   R := MultiplicatorRing(I);    // order for composite p   
                                 // if p is squarefree then this is
                                  // the correct result!
    
   O := Ep + R;
   
   return O;
end intrinsic;

intrinsic Optimize(p :: RngIntElt, O :: RngOrd) -> RngOrd
{computes a better representation of the maximal order}
   assert Discriminant(O) eq p^(p-2);

   SetOrderMaximal(O, true);

   test,Os := OptimizedRepresentation(O);  // takes a really long time

   return Os,test;
end intrinsic;

intrinsic GetUnits(p :: RngIntElt, O :: RngOrd) -> .
{computes the unit group of O (or a subgroup)}

   // SetKantVerbose("WILDAN", 3);
   C := ClassGroup(O:Bound := 1200);
   assert not IsDivisibleBy(#C,p);
   lp := Decomposition(O,p);
   U, m := SUnitGroup([lp[1][1],lp[2][1]]);
   
   return U,m;
end intrinsic;

intrinsic pAdicDefs (p :: RngIntElt) -> .
{defines the cyclotomic extension of Qp by pth roots of unity}
   eis := Evaluate(Qt!CyclotomicPolynomial(p),T+1);   // Eisenstein polynomial
   Qp := pAdicField(p: Precision := 50);
   Kp<t> := LocalField(Qp, eis, 20*(p-1));   // t = zeta_p - 1.	
	
   return t;
end intrinsic;

intrinsic pAdicRoot (p :: RngIntElt, r :: RngIntElt, t :: FldLocElt) -> FldLocElt
{finds a pth root of etap in Kp}
	etaval := &+[respow(p,i,p-r-1) : i in [1..(p-1)]];
	etap := &*[(1-(t+1)^i)^respow(p,i,p-r-1) : i in [1..p-1]]/t^etaval;
   	up := Root(etap,p : Precision := 15*(p-1));
   	wp1 := up*t^(etaval div p);

	return wp1;
end intrinsic;

intrinsic Embed (p :: RngIntElt, t :: FldLocElt, wp :: FldLocElt) -> .
{ returns the embedding of K into Kp and the autom. of Kp of order p-1}
	prim := PrimitiveRoot(p);
   	Kp := Parent(t);

   	function delp(i,x)
		e := Eltseq(x);
		v := Valuation(x);
  	
  		return &+[elt<Kp|e[j]*((t+1)^(respow(p,prim,i))-1)^(v+j-1)>: 
           j in [1..#e]];
   	end function;

   	function coerce(x)
		assert RelativePrecision(x) ge Precision(Kp);
	
		v := Valuation(x);
		e := Eltseq(x);
	
		return t^v * elt<Kp|e>;
   	end function;

   	function phi(x)
		e := Eltseq(x);
		valseq := [Valuation(e[j],p) : j in [1..#e] | not (e[j] eq 0)];
		if not (valseq eq []) then
			v := Min(valseq);
			prec := Precision(Kp)-v*(p-1);
			Lp := ChangePrecision(Kp,prec);
		end if;
	
		total := 0;
		for j in [1..#e] do
			if not (e[j] eq 0) then
				total +:= elt<Lp | e[j]*(Lp.1+1)^(j-1), prec>;
			end if;
		end for;
	
		if total eq 0 then
			return 0;
		else return coerce(total);
		end if;
   	end function;

   	function embed(i,x)
		e := Eltseq(x);
		K<z> := Parent(e[1]);
		taue := [ z^(i*(j-1))*e[j] : j in [1..#e]];
	
		return &+[phi(taue[j])*(wp)^(j-1) : j in [1..#e]];
   	end function; 
   	
   	return embed, delp;
end intrinsic;

intrinsic pAdicElt (p :: RngIntElt, wp :: FldLocElt, t :: FldLocElt,
   w :: RngUPolResElt, del) -> .
{modify wp by a root of unity in order to have the embedding commute with Galois}  

	embed,delp := Embed(p,t,wp);
   	while not (delp(1,wp) eq embed(0,del(w))) do
      	wp *:= t+1;
      	embed,delp := Embed(p,t,wp);
   	end while;
	
   	return wp;
end intrinsic;

intrinsic EmbedTrace (p :: RngIntElt, tr :: RngUPolResElt, t :: FldLocElt, wp :: FldLocElt) 
-> SeqEnum
{returns all the embeddings of tr in Kp}
   
   embed,_ := Embed(p,t,wp);
   v := [embed(i,tr) : i in [0..(p-1)]];

   return v;
end intrinsic;

intrinsic ExpandUnit (p :: RngIntElt, U, O :: RngOrd, g :: RngUPolElt, 
   a :: RngIntElt, m) -> RngIntElt
{returns a p-unit in K generating the unramified ideal above p in F}
   
   F<x> := NumberField(O);
   ug := F!O!(U.((p+3) div 2)@m);
   u := hom<F->F|F.1*a>(ug);
   
   return u;
end intrinsic; 


intrinsic Answer (p :: RngIntElt, u :: FldNumElt, v :: SeqEnum, t
:: FldLocElt) -> .
{returns N^1 of the p-unit u and checks to see if it is a pth power locally}
	
   useq := [hom<Parent(u)->Parent(t)|v[i]>(u) : i in [1..p]];
   ss := [elt<Parent(t)|ChangePrecision(useq[i],p+1)> : i in [1..p]];
   assert not ((&+[Valuation(useq[i]) : i in [1..p]] mod p) eq 0);   

   N1 := &*[ss[i+1]^i : i in [1..p-1]];
   assert (Valuation(N1) mod p) eq 0;
   N1t := N1/t^Valuation(N1);

   if N1t^2 eq 1+O(t^(p+1)) then ans := false;
   else ans := true;
   end if;
  
   return ans, ss, N1, N1t;
end intrinsic;