#_f2:= X^4 - +3*X^3-X^2 + 7;
#_k:= pAdicField(2,30);
#_g:= X^2 + 2*X + 2;
#_ek:= TotallyRamifiedExtension(_k,_g);
#_ekx:= PolynomialAlgebra(_ek);
# _f2:= Coerce(_ekx,_f2);
dd_realPrec:=20;




#_Z2q:= Quotient(pAdicRing(2,30),2^30);
#_Z2qx:= PolynomialAlgebra(_Z2q);
#_x:= _Z2qx.1;
#_f:= _x^6+6*_x^4+2^5*_x^3+12*_x^2-3*2^6*_x+33*2^3;

#_Z3:= pAdicRing(3,30);
#_Z3q:= Quotient(_Z3,3^30);
#_Z3qx:= PolynomialAlgebra(_Z3q);
#_y:= _Z3qx.1;
#_g:= _y^8+4*_y^6+6*_y^4+7*_y^2+9*_y+13;


#//freeze;

#/////////////////////////////////////////////////////////////////////////////
#// Automorphisms are in class_field.m now
#/////////////////////////////////////////////////////////////////////////////
#/*
#declare attributes Map:LocalAutoGenerators;
#declare attributes Map:LocalRingMap;

#function eval_automorphism(x, gens, pos)
#    if pos eq #gens then
#	return x;
#    end if;

#   R := Parent(x);
#    P := PolynomialAlgebra(R);
#    return Evaluate(
#	P ! [eval_automorphism(i, gens, pos + 1) : i in Eltseq(x)], gens[pos]
#    );
#end function;

#intrinsic Automorphisms(L::RngPad) -> SeqEnum
#{Return the automorphisms of the local ring L. Automorphisms are returned as a sequence of maps}
#    if L eq PrimeRing(L) then
#	curr_gens := [L ! 1];
#	curr := map<L -> L | x :-> eval_automorphism(x, curr_gens, 1)>;
#	curr`LocalAutoGenerators := curr_gens;
#	return [curr];
#    end if;

#    R := BaseRing(L);
#    maps := Automorphisms(R);
#    f := DefiningPolynomial(L);

#    res := [];
#    if IsEisenstein(f) then
#	P := PolynomialAlgebra(L);
#	for m in maps do
#	    gens := ChangeUniverse(m`LocalAutoGenerators, L);
#	    for r in Roots(P ! [m(x) : x in Eltseq(f)]) do
#		curr_gens := Insert(gens, 1, r[1]);
#		curr := map<L -> L | x :-> eval_automorphism(x, curr_gens, 1)>;
#		curr`LocalAutoGenerators := curr_gens;
#		Append(~res, curr);
#	    end for;
#	end for;
#    elif IsInertial(f) then
#	for m in maps do
#	    gens := ChangeUniverse(m`LocalAutoGenerators, L);
#	    r := L.1;
#	    for i in [1..Degree(L)] do
#		curr_gens := Insert(gens, 1, r);
#		curr := map<L -> L | x :-> eval_automorphism(x, curr_gens, 1)>;
#		curr`LocalAutoGenerators := curr_gens;
#		Append(~res, curr);
#		r := GaloisImage(r, 1);
#	    end for;
#	end for;
#    else
#	error "Error: Unsupported defining polynomial";
#    end if;

#    return res;
#end intrinsic;

#function eval_field_automorphism(x, m)
#    F := Parent(x);
#    R := RingOfIntegers(F);
#    pi := UniformizingElement(F);

#    xv := Valuation(x);
#    xu := ShiftValuation(x, -xv);
#    return (F ! m(UniformizingElement(R)))^xv * F ! m(R ! xu);
#end function;

#function construct_field_map(L, m)
#    res := map<L -> L | x :-> eval_field_automorphism(x, m)>;
#    res`LocalRingMap := m;
#    return res;
#end function;

#intrinsic Automorphisms(L::FldPad) -> SeqEnum
#{Return the automorphisms of the local ring L. Automorphisms are returned as a sequence of maps}
#    maps := Automorphisms(RingOfIntegers(L));
#    return [construct_field_map(L, m) : m in maps];
#end intrinsic;

#function auto_eq(m1, m2)
#    return m1`LocalAutoGenerators eq m2`LocalAutoGenerators;
#end function;

#function auto_mult(m1, m2)
#    gens1 := m1`LocalAutoGenerators;
#    gens2 := m2`LocalAutoGenerators;
#    gens := [];
#    for i in [1..#gens1] do
#	Append(~gens, eval_automorphism(gens1[i], gens2, 1));
#    end for;

#    L := Universe(gens1);
#    m := map<L -> L | x :-> eval_automorphism(x, gens, 1)>;
#    m`LocalAutoGenerators := gens;
#    return m;
#end function;

#intrinsic AutomorphismGroup(L::RngPad) -> GrpPerm
#{Return the automorphism group of L as a permutation group}
#    maps := Automorphisms(L);

#    G, m1 := GenericGroup(
#	maps : Mult := auto_mult, Eq := auto_eq, Id := maps[1]
#    );
#    m2, G2 := CosetAction(G, sub<G|>);
#    m := Inverse(m2) * m1;
#    return
#	G2,
#	map<
#	    G2 -> PowerStructure(Map) | g :-> m(g), h :-> h @@ m
#	>;
#end intrinsic;

#intrinsic AutomorphismGroup(L::FldPad) -> GrpPerm
#{Return the automorphism group of L as a permutation group}
#    G, m := AutomorphismGroup(RingOfIntegers(L));

#    return
#	G,
#	map<
#	    G -> PowerStructure(Map) |
#	    g :-> construct_field_map(L, m(g)),
#	    h :-> h`LocalRingMap @@ m
#	>;
#end intrinsic;
#
#///////////////////////////////////////////////////////////////////////
#// general functions
#//////////////////////////////////////////////////////////////////////

#intrinsic BaseRing(L::FldPad) -> FldPad
#{The base field of L}
#  return BaseField(L);
#end intrinsic;


#//---------------------------------------------------------------------
#// FACTORIZATION
#//---------------------------------------------------------------------

# default precision
_realPrec := 10;

#///////////////////////////////////////////////////////////////////////
#// some functions for polynomials that probably should be somewhere
#// else

_ChangePrecisionCoeff:= function(PR,n)
  local Co,Co2;
  Co:= CoefficientRing(PR);
  Co2:= ChangePrecision(Co,n);
  return PolynomialAlgebra(Co2);
end;




IsMonic:= function(f)
    return LeadingCoefficient(f)= One(CoefficientRing(Parent(f)));
end;

#Computes a polynomial g such that g mod Lm[i] = La[i], Lm and La are lists of Polynomials
ChineseRemainder:= function(Lm,La)
local PR,Mod,CR,i,gcd,e2,e1;

PR := Parent(Lm[1]);
Mod := Lm[1];
CR := La[1];
for i in [2..Length(Lm)] do
    gcd := XGCD(Mod,Lm[i]);

    if Base(gcd)<> One(PR) then
       return Error( "Error, ChineseRemainderTheorem: modules must be coprime.");
    fi;

    e2 := Mod*(gcd.ext1);
    e1 := Lm[i]*(gcd.ext2);

    Mod := Mod*Lm[i];
    CR := (e1*CR+e2*La[i]) mod Mod;
  od;
  return CR;
end;

PGcd:=function(f,g)
  local PR,R,pi,prec,S,H,L,h,loss,l,i,H,p,k;
#Print("PGcd g: ",g,"\nf ",f,"\n");
#// only works for quotient rings.
  if Degree(f) < 0 or Degree(g) < 0 then return [One(Parent(f)),0]; fi;
#  ////////////////assert IsMonic(f) or IsMonic(g);
  #Print("f: ",f);
  #Print("\ng: ",g,"\n");
  PR := Parent(f);
  R:= CoefficientRing(PR);
  pi:= UniformizingElement(R);
  prec := Precision(R);
  S := SylvesterMatrix(f,g);
  H := EchelonForm(S);
#Print("\nPGcd HNF",H,"\n");
  H:=H.base;
  L := ElementToSequence(H[Rank(H)]);
#Print("L ",L);
  h := Coerce(PR,Reverse(L));
  l := ElementToSequence(h);
  k:= [];
  for i in [1..Length(l)] do
      k[i]:= Valuation(l[i]);
  od;
  loss:=Minimum(k);
  H:= ElementToSequence(h);
  p:= pi^loss;
  for i in [1..Length(H)] do
      H[i]:= Div(H[i],p);
  od;
  h := Coerce(PR,H);
#Print("\nPGcd ",h," loss ",loss,"\n");  
  return [h, loss];
end;


#///////////////////////////////////////////////////////

is_squarefree:= function(f)
   local val,df,L,gcd,loss;
  if Degree(f) <= 0 then
    return [true, f];
  fi;

 # vprint RoundFour,4:" -- disc";
 # vtime RoundFour,4:
 if Is(Type(f),elt-alg^pol/res^pad) or Is(Type(f),elt-alg^pol/ord^pad) then 
 val :=  Valuation(Discriminant(f));
  if val < Precision(CoefficientRing(Parent(f))) then
    #vprint RoundFour,3:" -- polynomial is squarefree by discriminant";
   return [true, f, 0];
  fi;
  fi;
  df := Derivative(f);
  #vprint  RoundFour,4: " -- gcd";
  #vtime RoundFour,4:
  L := PGcd(f,df);
  gcd:= L[1];
  loss:= L[2];
  if Degree(gcd) < 1 then
    #vprint RoundFour,3:" -- polynomial is squarefree by gcd";
    return [true, f, loss];
  else
    return [false,Div(f,gcd), loss];
  fi;

end;

#//TODO ???
#//intrinsic IsSquareFree(f::RngUPolElt[RngPadResExtElt])
#//{return
#//  is_squarefree(f)



#///////////////////////////////////////////////////////
#// functions for handling polynomials with denominators

#/////////////////////////////////////////////////////////////
#// in the round four elements are represented by polynomials
#// and the exponent of pi in the denomiantor
#// we also store additional information in the structure
#theta has to be a record with .elt and .denval components

# 'cancels out' powers of the uniformizer pi in the polynomial
# and its denominator.

elt_reduce := function(theta)
local seq,vallist,i,minval,pi;
  if theta.elt = Coerce(Parent(theta.elt),0) then
    theta.denval:= 0;
    return theta;
  elif theta.elt<>"null" then
    seq := ElementToSequence(theta.elt);
    vallist:= [];
    for i in [1..Length(seq)] do
        vallist[i]:= Valuation(seq[i]);
    od;
    minval := Minimum(vallist);
    if minval > 0 and theta.denval>= minval then
      pi := UniformizingElement(CoefficientRing(Parent(theta.elt)));
      theta.elt := Div(theta.elt,pi^minval);
      theta.denval := theta.denval - minval;
    fi;
  fi;
return theta;
end;

# constructor
eltRec := function(phi, denval)
local theta;
 #   ELT :=    recformat
   #           <elt:RngUPolElt,       // the numerator of the element
  #             denval:IntegerRing(),   // the denominator of elt is pi^den,
   #            chi:RngUPolElt,       // its characteristic polynomial
     #          passes_hensel_test,   // true if it passes the hensel test
    #           nu:RngUPolElt,        // irreducible factor of chi mod pi
      #         f:IntegerRing(),         // deg(nu)
       #        res_fact,             // factorization of chi mod pi
        #       passes_newton_test,   // true if it passes the newton test
         #      v_star,               // v_star valuations of elt
          #     e:IntegerRing()          // numerator of vstar
           #   >;
    theta := rec(elt := phi, denval:= denval, chi:="null", passes_hensel_test:="null", nu:="null",f:="null", res_fact:="null",passes_newton_test:="null",v_star:="null",e:="null");
    theta :=elt_reduce(theta);
    return theta;

end;

#//////////////////////////////////////////////////
#// Find the inverse of some (invertible) element theta modulo Phi.
#// the inverse is given as res/pi^denval
elt_inverse := function(theta, Phi)

local pi,P,R,M,L,E,T,res2,gcd,denval,uni;
    pi := UniformizingElement(CoefficientRing(Parent(Phi)));
    P := Parent(theta.elt);
    R := CoefficientRing(P);
    M := SylvesterMatrix(theta.elt, Phi);
    L:=EchelonForm(M);
    E:= L.base;
    T:= L.ext1;
    res2 := (Coerce(P,Reverse(ElementToSequence(ExtractBlock(T, NumberOfRows(T), 1,1,Degree(Phi)))))) mod Phi;
    gcd := Coerce(CoefficientRing(E),E[NumberOfRows(E)][NumberOfColumns(E)]);
    denval := Valuation(gcd);
    uni := Div(gcd,pi^denval);
    res2 := Div(res2,uni);
    denval := denval - theta.denval;
    if denval < 0 then
      res2 := res2*pi^(-denval);
      denval := 0;
    fi;
    return eltRec(res2,denval);
end;

#//////////////////////////////////////////////////
# compute the v_star valuation from a given characteristic polynomial f
v_star := function(f)
local coeffs,L,i, res2;
    if IsMonic(f)= false then
       Print("f not monic\n");
       return;
    fi;
    coeffs := Reverse(ElementToSequence(f));
    Butfirst_(coeffs);
    L:=[];
    for i in [1..Length(coeffs)] do
        L[i]:= (Valuation(coeffs[i]))/i;
    od;
    res2 := Minimum(L);
    return res2;
end;


#////////////////////////////////////////////
#// computation of characteristic polynomials

#////////////////////////////////////////////////////////////////
#// compute the characteristic polynomials using newton relations
#// (cohen 1, p 161)
#// the traces_from_poly of tyhe reference polynomials Phi are
#// precomputed and applied in th computation of char polys wrt
#// Phi.

traces_from_poly := function(Phi)
local n,R,L,inc,tmp_R,tmp_PR,tmp_Phi,t,S,k,Sk,i,char;
  n := Degree(Phi);
  R := CoefficientRing(Parent(Phi));
  if Is(Type(R),ord^ser) or Is(Type(R),res^ser) then
     char := Coerce(R,Characteristic(CoefficientRing(R)));
  else char:= -1;
  fi;
  #// precision loss occurs in poly_from_traces by division
  if n > 0 then
    L:=[];
    for i in [1..n] do
        if Coerce(R,i) = char then
           L[i]:= 0;
        else L[i]:= Valuation(Coerce(R,i));
        fi;
    od;
    inc := Maximum(L);
    inc := inc^2;
  else
    inc := 0;
  fi;
  if not Precision(R)=INFTY then
  tmp_R := ChangePrecision(R,Precision(R)+inc);
  else tmp_R:= R;
  fi;
  tmp_PR := PolynomialAlgebra(tmp_R);
  tmp_Phi := Coerce(tmp_PR,Phi);

  t := function(i)
    if i < 0 then
      return 0;
    else
      return Coefficient(tmp_Phi,i);
    fi;
  end;

  S := [];

  for k in [1..n] do
    Sk := -k*t(n-k);
    for i in [1..Minimum(n,k-1)] do
      Sk := Sk - t(n-i)*S[k-i];
    od;
    Append_(S,[Sk]);
  od;
  return S;
end;


poly_from_traces := function(S)
local n,Tn,k,Tnk,i;
  n := Length(S);
  Tn := [];

  for k in [1..n] do
    Tnk := -S[k];
    for i in[1..k-1] do
      Tnk := Tnk - Tn[i]*S[k-i];
    od;
    Tnk := Div(Tnk,k);
    Append_(Tn,[Tnk]);
  od;
  Tn := Reverse(Tn);
  Append_(Tn,[Coerce(Parent(S[1]),1)]);
  return Coerce(PolynomialAlgebra(Parent(S[1])),Tn);
end;


chi := function(Phi, traces_Phi, theta)
local ret,k, L,tmp_PR,tmp_theta,tmp_Phi,n,S,pow2,i,Si,j,M,S,li,T,nu;

  tmp_PR := PolynomialAlgebra(Parent(traces_Phi[1]));
  #L:= ElementToSequence(theta);
  #Print(theta,"\n");
  #Print(L,"\n",Type(L[2]),"\n");
 #M:= [];
     # for i in [1..Length(L)] do
	 #    nu:= ElementToSequence(L[i]);
             
	  #   if Length(nu)>1 then
	  #      li:= false;
	   #     for j in [2..Length(nu)] do
		#    if nu[j] <>0 then
	#	       li:= true;
		#    fi;
		# od;
	      #  if li then
	      #     T:= [];
	      #     for j in [1..Length(nu)] do
	     #         T[j]:=Coerce(IntegerRing(),nu[j]);
	    #       od;
	      #     M[i]:= Element(CoefficientRing(tmp_PR),T);
	     #   else
              #     M[i]:=Coerce(IntegerRing(),L[i]);
		#fi;
	    # else
	   #  	     M[i]:=Coerce(IntegerRing(),L[i]);
	#  fi;
      #od;

  tmp_theta := Coerce(tmp_PR,theta);
  tmp_Phi := Coerce(tmp_PR,Phi);
#Print("tmp_PR: ",tmp_PR,"\nCoefficientRing(tmp_PR): ",CoefficientRing(tmp_PR),"\nType(theta): ",Type(theta),"\nCoefficientRing(Parent(theta)): ",CoefficientRing(Parent(theta)),"\n");
  n := Degree(Phi);
  S := [];
  pow2 := 1;
  for i in [1..n] do
    pow2 := pow2*tmp_theta mod tmp_Phi;

    Si := 0;
    Si := n*Coefficient(pow2,0);
    for j in [1..Degree(pow2)] do
      if Coefficient(pow2,j) <> 0 then
        Si :=Si+traces_Phi[j]*Coefficient(pow2,j);
      fi;
    od;
    Append_(S,[Si]);
  od;
  ret :=  Coerce(Parent(Phi),poly_from_traces(S));
  return ret;
end;

pow_chi:= function(Phi,theta)
local A2,chi,L,j,A,chir,a,b, elseq,pi,poly,prec,n,i,j,E,ech,phi,pa,T,nu,ra,k,loss,p,inv,denv,mult;
phi:= theta.elt;
denv:= theta.denval;
pa:= Parent(Phi);
prec:= Precision(CoefficientRing(pa));
pi:= UniformizingElement(CoefficientRing(pa));
n:= Degree(Phi);
A:= [];
for i in [1..n+1] do
    a:= pi^(denv*(n-i+1));
    b:= phi^(i-1);
    a:= Coerce(Parent(b),a);
    elseq:= ElementToSequence(a*b mod Phi);
    for  j in [Length(elseq)+1..n+1] do
         Append_(elseq,[Coerce(CoefficientRing(pa),0)]);
    od;
    A:=Append(A,elseq);
od;
for i in [1..n+1] do
    L:=[];
    for j in [1..n+1] do
        if i = j then 
           L[j]:= pi^prec;
        else L[j]:= Coerce(pa,0);
        fi;
    od;
    A:=Append(A,L);
od;
A:= Matrix(CoefficientRing(pa),2*(n+1),n+1,A);
ech:= EchelonForm(A);
E:= ech.base;
T:= ech.ext1;
ra:= Rank(E);
for i in [ra..2*(n+1)] do
    nu:= true;
    for j in [1..Length(E[i])] do
        if E[i][j]<> Coerce(Parent(E[i][j]),0) then
           nu:=false;
           j:= Length(E[i]);
        fi; 
   od;
if nu then chi:= Coerce(Parent(Phi),List(T[i]));
   if chi <> Element(Parent(Phi),0) then
       chir:=chi;
   fi;
fi;
od;
elseq:= ElementToSequence(chir);
k:= [];
poly:=[];
  for i in [1..Length(elseq)] do
      poly[i]:= elseq[i];
      k[i]:= Valuation(poly[i]);
  od;
loss:=Minimum(k);
p:= pi^loss;
chir:= Element(Parent(Phi),poly);
chir:= Div(chir,p);
inv:= IsInvertible(LeadingCoefficient(chir));
if inv.base then
   inv:= inv.ext1;
else return Error("pow_chi not monic");
fi;
chir:= chir*Coerce(pa,inv);
return [chir,loss];
end;    


#{Find char poly of theta * pi^-theta_denval}
elt_chi := function(Phi, traces_Phi, theta)
local PR,R,tmp_R,tmp_PR,res2,pi;
     if Is(Type(Phi),elt-alg^pol/ord^ser) or Is(Type(Phi),elt-alg^pol/res^ser) then
       return  pow_chi(Phi,theta)[1];
    fi;
    #vprint RoundFour,5:" -- char poly";
    #vprint RoundFour,6:" -- of",theta`elt;
    #vprint RoundFour,6:" -- with denominator pi ^",theta`denval;
    PR := Parent(Phi);
    R := CoefficientRing(Parent(Phi));
    if not Precision(R)=INFTY then
    tmp_R := ChangePrecision(R,Precision(R)+4*theta.denval*Degree(Phi));
    else tmp_R:= R;
    fi;
    tmp_PR := PolynomialAlgebra(tmp_R);
    #vtime RoundFour,5:
    #Print("Phi: ",Phi,"\nCoefficientRing(Parent(Phi)):",CoefficientRing(Parent(Phi)),"\n");
    #Print("traces_Phi: ",traces_Phi,"\nType(traces_Phi[1]): ",Type(traces_Phi[1]),"\n");
    #Print("theta.elt:",theta.elt,"\nCoefficientRing(Parent(theta.elt)):",CoefficientRing(Parent(theta.elt)),"\n\n");
    res2 := chi(Phi, traces_Phi, theta.elt);
    #//"vstar",v_star(res);
    #//"denval",theta`denval;
    if v_star(res2) < theta.denval then
      return Error("Factorization:Insufficient precision in characteristic polynomial computation");
    fi;
    res2 := Coerce(tmp_PR,res2);
    pi := UniformizingElement(tmp_R);
    res2 := Evaluate(res2, pi^theta.denval * tmp_PR.1);
    res2:= Div(res2,LeadingCoefficient(res2));
    res2 := Coerce(PR,res2);
     if IsMonic(res2)= false then
       return Error("f not monic at elt_chi\n");
       return;
    fi;
#Print(" testing chi(theta)", elt_evaluate(Phi,eltRec(Coerce(PR,res2),0),theta),"\n");
    return Coerce(PR,res2);
end;


#// compute invariants of elt, war eigentlich procedure, die phi veraenderte
elt_invariants:= function(Phi, traces_Phi, phi)
local res_chi,fact,K,pc;

  if phi.chi="null" then
    if Is(Type(Phi),elt-alg^pol/ord^ser) or Is(Type(Phi),elt-alg^pol/res^ser) then
       pc:= pow_chi(Phi,phi);
       phi.chi:= pc[1];
       phi.chi:= Coerce(PolynomialAlgebra(ChangePrecision(CoefficientRing(Parent(pc[1])),Precision(CoefficientRing(Parent(pc[1])))+pc[2])),phi.chi);
    else
    phi.chi := elt_chi(Phi, traces_Phi, phi);
    fi;
  fi;
  phi.v_star := v_star(phi.chi);
  #Print("v_star:",phi.v_star,"    chi:",phi.chi,"\n");
  #Print("CoefficientRing(Parent(chi)): ", CoefficientRing(Parent(phi.chi)),"\n");
  phi.passes_newton_test :=
    phi.v_star = Valuation(TrailingCoefficient(phi.chi)) / Degree(phi.chi);
  phi.e := Denominator(phi.v_star);
  if phi.passes_newton_test or phi.v_star = 0 then
    if phi.res_fact = "null" then
      if IsField(CoefficientRing(Parent(phi.chi))) then
         K:=ResidueClassField(IntegerRing(CoefficientRing(Parent(phi.chi))));
      else
          K := ResidueClassField(CoefficientRing(Parent(phi.chi)));
      fi;

      #Print("K: ",K,"phi.chi: ", phi.chi,"\n");
      res_chi := Coerce(PolynomialAlgebra(K),phi.chi);
     # Print("Parent(res_chi): ",Parent(res_chi),"res_chi: ",res_chi,"\n");
      fact := Factorization(res_chi);
      #Print("fact: ",fact," Parent(fact[1][1]): ",Parent(fact[1][1]),"\n");
      phi.res_fact := fact;
    fi;
    phi.passes_hensel_test := (Length(phi.res_fact) = 1);
    if phi.passes_hensel_test then
      phi.nu := Coerce(Parent(phi.chi),phi.res_fact[1][1]);
      phi.f := Degree(phi.nu);
   fi;
  else
    phi.passes_hensel_test := true;
  fi;
 return phi;
end;

#// ??? !!! TODO ft(k2,L[14][3]);

#// returns theta`elt(tau`elt/pi^tau`denval)/pi^theta`denval mod Phi
elt_evaluate := function(Phi, theta, tau)
local PR,R,inc,tmp_R,tmp_PR,pi,tmp_Phi,tmp_theta,tmp_tau,res2,L,se,i,res_denval;
  PR := Parent(Phi);
  R := CoefficientRing(PR);
  inc := Degree(theta.elt)*tau.denval;
  if not Precision(R)=INFTY then
  tmp_R := ChangePrecision(R,Precision(R)+inc);
  else tmp_R:= R;
  fi;
  tmp_PR := PolynomialAlgebra(tmp_R);
  pi := UniformizingElement(tmp_R);
  tmp_Phi := Coerce(tmp_PR,Phi);
  tmp_theta :=Coerce(tmp_PR,theta.elt);
  tmp_tau := Coerce(tmp_PR,tau.elt);


  res2 := tmp_theta*pi^inc;
  L:= [];
  se:= ElementToSequence(res2);
  for i in [1..Degree(res2)+1] do
        L[i]:= Div(se[i],pi^(tau.denval*(i-1)));
  od;
  res2 := Coerce(tmp_PR,L);
  res2 := Evaluate(res2, tmp_tau) mod tmp_Phi;
  res_denval := theta.denval+inc;
  return eltRec(Coerce(PR,res2), res_denval);
end;

#// elt_evaluate
#// and transfer invariants
#// useful when changing reference polynomial
elt_evaluate_trans := function(Phi, traces_Phi, theta, tau)
local res2;

  res2 := elt_evaluate(Phi, theta, tau);
  if theta.chi <> "null" then
    res2.chi := theta.chi;
  fi;
  if theta.res_fact <>"null" then
    res2.res_fact := theta.res_fact;
  fi;
  res2:=elt_invariants(Phi, traces_Phi, res2);
  return res2;
end;



#//////////////////////////////////////////////////////////////////////////
#//
#// returns a list of exponents res such that
#// v*(pi^res[1]*pool[1]^res[2]*...*pool[m]^res[m+1]) eq val
#// where res[i]<ceils[i-1] (for i>=2) and
#// Denominator(val) | E
calculate_exponents := function(pool, ceils, val, E)
local res2,rem_val,len,now,d,tmpval,tmpool,mud,rcr,sol;
  res2 := [];
  rem_val := val;
  len := Length(pool);
  now := len;

  d := E;

  while now > 0 do

    d := d/(ceils[now]);
    tmpval := rem_val * d;
    tmpool := pool[now] *d;
    mud := Denominator(tmpool);
    rcr := ResidueClassRing(Coerce(IntegerRing(),mud));
    tmpval := tmpval * mud;
    tmpool := tmpool * mud;
    sol := Coerce(IntegerRing(),((Coerce(rcr,tmpval))/(Coerce(rcr,tmpool))));

    rem_val := rem_val - sol * pool[now];
    Append_(res2, [sol]);

    now := now - 1;

  od;

  Reverse_(res2);
  Insert_(res2,rem_val,1);
  return res2;
end;

#//{Returns an element poly, such that
#// poly has v_star valuation vstar}
elt_with_v_star := function(Phi, P, P_v_star, S, E, vstar)
local exps,poly_val,poly,i,pi;

	exps := calculate_exponents(P_v_star, S, vstar, E);
	poly_val := Coerce(IntegerRing(),exps[1]);
	poly := Coerce(Parent(Phi),1);
	if Length(P) <> 0 then
	    for i in [1..Length(P)] do
		 poly:= poly*P[i] ^ exps[i + 1];
	    od;
	fi;


        if poly_val > 0 then
          pi := UniformizingElement(CoefficientRing(Parent(Phi)));
          poly := poly*pi^poly_val;
          poly_val := 0;
        fi;
        return eltRec(poly, -poly_val);
end;

#//////////////////////////////////////////////////////////////////////////
#// the denominator of theta is pi^theta`denval
#// if theta`chi is not assigned or theta`denval gt 0 then chi_theta is computed

hensel_lift_factorization_sub := function(Phi, theta)
local PR,R,tmp_R,tmp_PR,K,P,pi,tmp_Phi,traces_tmp,Phi,tmp_theta,facts,i,j,lifted_facts,lf,                 eval_facts,res2,traces_tmp_Phi,elm,L,elmi;
    PR := Parent(Phi);
    R := CoefficientRing(Parent(Phi));
    if not Precision(R) =INFTY then
    tmp_R := ChangePrecision(R,Precision(R)+Degree(Phi)*(theta.denval));
    else tmp_R:= R;
    fi;
    tmp_PR := PolynomialAlgebra(tmp_R);
    K := ResidueClassField(tmp_R);
    P := PolynomialAlgebra(K);
    pi := UniformizingElement(tmp_R);
    tmp_Phi := Coerce(tmp_PR,Phi);
    traces_tmp_Phi := traces_from_poly(tmp_Phi);
    tmp_theta := theta;
    if theta.denval <> 0 or theta.chi ="null" then
      tmp_theta.chi := elt_chi(tmp_Phi, traces_tmp_Phi, tmp_theta);
    fi;
#Print(" testing chi(theta)", elt_evaluate(tmp_Phi,eltRec(tmp_theta.chi,0),tmp_theta),"\n");
    if theta.res_fact <>"null" then
      tmp_theta.res_fact := theta.res_fact;
    fi;
    tmp_theta := elt_invariants(tmp_Phi, traces_tmp_Phi,tmp_theta);
    facts:=[];
    if tmp_theta.res_fact<>"null" then
        for i in [1..Length(tmp_theta.res_fact)] do
     facts[i]:=Coerce(tmp_PR,Coerce(P,(tmp_theta.res_fact[i][1]^tmp_theta.res_fact[i][2])));
        od;
    fi;
    if Length(facts)=1 then
      if facts[1] = tmp_theta.chi then
         L:=[];
         L[1]:= [tmp_theta.chi,0];
         return L;
      else
      facts:=Append(facts,[Coerce(Parent(facts[1]),1)]);
      fi;
    fi;
    #Print(facts);
    lifted_facts := HenselLift(tmp_theta.chi, facts);
    lf:=lifted_facts[1];
    for i in [2..Length(lifted_facts)] do
        lf:= lf*lifted_facts[i];
    od;
    if Is(Type(R),ord^ser) or Is(Type(R),res^ser) then
       # assuming precision 20
       elm:=(ElementToSequence(lf-tmp_theta.chi));
        if elm <> [] then
        for i in [1..Length(elm)] do
            elmi:=ElementToSequence(elm[i]);
                for j in [1.._realPrec] do
                   if elmi[j]<>Element(Parent(elmi[i]),0) then return Error("assertion violated\n");fi;
                od;
        od;
        fi;
    else if (lf- tmp_theta.chi)<>0 then return Error("assertion violated\n"); fi;
    fi;

    eval_facts:= [1..Length(lifted_facts)];
#Print("henseltest lifted_facts ",lifted_facts);
    Apply_(eval_facts,x->elt_evaluate(tmp_Phi,eltRec(lifted_facts[x],0),tmp_theta));
#Print("henseltest eval_facts ",eval_facts,"\n");
    #vprint  RoundFour,4: " -- gcd";
    #vtime RoundFour,4:
    res2:=[1..Length(eval_facts)];
    Apply_(res2,x->PGcd(Phi,Coerce(PR,(eval_facts[x].elt))));
    return res2;
end;


#// the denominator of theta is pi^theta_denval
#// if chi_theta eq 0 or theta_denval gt 0 then chi_theta is computed
hensel_lift_factorization_chi2 := function(Phi, theta)
local flag,PR,R,res2,inc,resProd,i,tmp_R,tmp_theta,res2,tmp_PR,re;
    #Print("hensel_lift_factorization_chi2");
    flag := 1;

    PR := Parent(Phi);
    R := CoefficientRing(Parent(Phi));
   # vprint RoundFour,5:" == Hensel lifting to precision", Precision(R),"...";
    res2 := hensel_lift_factorization_sub (Phi, theta);
    #// if the precion was not sufficient, repeat this step with a higher
    #// precision.
    inc := 1;
    resProd:= res2[1][1];
       while (IsList(resProd)) do
           resProd:= resProd[1];
       od;
    for i in[2..Length(res2)] do
        if IsList(res2[i]) then
        re:=res2[i][1];
	while IsList(re) do
	    re:= re[1];
	od;
        resProd:= resProd*re;
        fi;
    od;
#Print("\nResProd ",resProd,"\n");
    while (flag = 1 and Element(Parent(Phi),0) <>((resProd)-Phi)) do
      #Print("difference: ",resProd-Phi,"\n");
      #//print inc;
      if not Precision(R) =INFTY then
      tmp_R := ChangePrecision(R,Precision(R)+4*Degree(Phi)*inc);
      else tmp_R:= R;
      fi;
      #vprint RoundFour,5:" == ... trying Hensel lift again, precision increased by",
       #                  4*Degree(Phi)*inc,"...";
      tmp_PR := PolynomialAlgebra(tmp_R);
      tmp_theta := eltRec(Coerce(tmp_PR,theta.elt), theta.denval);

      if theta.res_fact<>"null" then

        tmp_theta.res_fact := theta.res_fact;
      fi;
      #Print("henselLift: prec ",Precision(tmp_R));
      res2:= hensel_lift_factorization_sub(Coerce(tmp_PR,Phi),tmp_theta);
      for i in [1..Length(res2)] do
         res2[i][1] := [Coerce(PR,res2[i][1])];
      od;
      resProd:= res2[1][1];
      while (IsList(resProd)) do
           resProd:= resProd[1];
       od;
      for i in[2..Length(res2)] do
         re:=res2[i][1];
	 while(IsList(re)) do
	    re:= re[1];
	 od;
         resProd:= resProd*re;
      od;
      inc :=inc+ 1;
    od;
    #vprint RoundFour,5:" == ... done";
    return res2;
end;


newton_lift_factorization_chi2 := function(Phi, theta)
local rs,gamma;
# //print "newton theta",theta,theta_denval;
#//"theta",theta;
  if theta.v_star = 0 and theta.passes_hensel_test=FALSE then
    #//print "newton_lift_factorization_chi2 hensel";
    return hensel_lift_factorization_chi2(Phi, theta);
  fi;

  rs := theta.v_star;
  #vprint RoundFour,5:" == flattening Newton polygon slope from",rs,"to 0";
  gamma := eltRec(Coerce(Parent(Phi),(theta.elt^Denominator(rs)) mod Phi),
               Numerator(rs)+theta.denval*Denominator(rs));
  return hensel_lift_factorization_chi2(Phi, gamma);

end;

#/////////////////////////////////////////////////////////////////////////////////

#//{Find a polynomial delta with Degree(delta)<Degree(nu_phi)
#// and v_star(delta-gamma)>0}
#// returns:
#// true, delta, xxx
#//   if such an element delta has been found
#// false, xxx, rho
#//   if during the search for delta an element that fails the
#//   hensel test has occured.

res_field_rep := function (Phi, traces_Phi, nu, gamma)
local pi,K,P,nur,F,PF,facts,i,a,b,delta,rho,Co,p,prec,g,L,M,S,T,j,U,li;

   if Degree(gamma.nu) = 1 then
     return [true, -TrailingCoefficient(gamma.nu), 0];
   fi;

   if gamma.denval = 0 then
     return [true, gamma mod nu, 0];
   fi;
   Co:= CoefficientRing(Parent(Phi));

  # if Is(Type(Co),res^pad) then
   #   p:= Size(ResidueClassField(Co).base);
    #  prec:= Coerce(IntegerRing(),(Precision(Co)/RamificationDegree(Co)));
    #  g:= DefiningPolynomial(Co);
    #  Co:= pAdicRing(p,prec);
    #  Co:= TotallyRamifiedExtension(Co,g);
    #  L:= ElementToSequence(g);

    #  M:= [];
    #  for i in [1..Length(L)] do
    #	     S:= ElementToSequence(L[i]);
#	     if Length(S)>1 then
#	        li:= false;
#	        for i in [2..Length(S)] do
#		    if S[i] <>0 then
#		       li:= true;
#		    fi;
#		 od;
#	        if li then
#	           T:= [];
#	           for j in [1..Length(S)] do
#	              T[i]:=Coerce(IntegerRing(),S[i]);
#	           od;
#	           M[i]:= Element(Co,T);
#	        else
#                   M[i]:=Coerce(IntegerRing(),L[i]);
#		fi;
#	     else M[i]:=Coerce(IntegerRing(),L[i]);
#	  fi;
#     od;
#      gamma.elt:= Coerce(PolynomialAlgebra(Co),M);
#   fi;
   pi := UniformizingElement(Co);
   K := ResidueClassField(Co);
   P := PolynomialAlgebra(K);
   nur := Coerce(P,nu);
   F := Extension(K,nur);
   PF:= PolynomialAlgebra(F);
   facts := Roots(Coerce(PF,(Coerce(P,gamma.nu))));
   for i  in [1..Length(facts)] do
     a := ElementToSequence(facts[i][1]);
     b := Coerce(PF,a);
     delta := Coerce(PolynomialAlgebra(Co),b);
     rho := eltRec(gamma.elt-Coerce(Parent(gamma.elt),pi)^gamma.denval*delta,gamma.denval);
     rho:=elt_invariants(Phi, traces_Phi, rho);
     if rho.passes_hensel_test=FALSE then
       return [false, 0, rho];
     fi;
     if rho.v_star > 0 then
       return [true, Coerce(Parent(Phi),delta), 0];
     fi;
   od;
   return Error ("Factorization: Round Four Panic, no representing element found.");
end;

#// x gives a residue class field extension of degree F_x
#// gamma/pi^gamma_denval gives a residue class field extension of
#// degree F_gamma
#// find a poltynomial that gives a residue class field extension of
#// degree LCM(F_x,F_gamma)
#  //print "kappa",kappa,"/",pi^kappa_denval;

ff_composition_gen := function(Phi, x, F_x, gamma)
local kappa,PR,R,denval,tmp_R,tmp_PR,x,pi,tmp_Phi,traces_tmp_Phi,tmp_gamma,K,F,pow,pow2,coeff,li,y,i;
  if gamma.f = LCM(F_x,gamma.f) then
    kappa := gamma;
    return kappa;
  fi;

  #//print "gamma",gamma;
  PR := Parent(Phi);
  R := CoefficientRing(Parent(Phi));
  denval := gamma.denval*gamma.f;
  if not Precision(R) = INFTY then
  tmp_R := ChangePrecision(R,Precision(R)+denval);
  else tmp_R:= R;
  fi;
  #//tmp_R := pAdicQuotientRing(Prime(R),Precision(R)+denval);
  tmp_PR := PolynomialAlgebra(tmp_R);
  y:=tmp_PR.1;
  pi := UniformizingElement(tmp_R);
  tmp_Phi := Coerce(tmp_PR,Phi);
  traces_tmp_Phi := traces_from_poly(tmp_Phi);
  tmp_gamma := Coerce(tmp_PR,gamma.elt);

  K := ResidueClassField(tmp_R);

  F := LCM(F_x,gamma.f);
  pow:=[1..F_x];
  Apply_(pow,x->(pi^(denval)*y^(x-1)) mod (pi^denval*tmp_Phi));
  pow2:=[1..gamma.f];
  Apply_(pow2,x->pi^(denval-(gamma.denval*(x-1)))*tmp_gamma^(x-1)mod (pi^denval*tmp_Phi));
  pow := Concatenation(pow,pow2);

  repeat
    coeff:=[1..Length(pow)];
    Apply_(coeff,x->Coerce(tmp_R,Random(K)));
    #//"coeff",coeff;
    li:= [1..Length(pow)];
    Apply_(li,x->coeff[x]*pow[x]);
    kappa:= li[1];
    for i in [2..Length(li)] do
       kappa := kappa+li[i];
    od;
    kappa:= eltRec(kappa mod tmp_Phi,denval);
    kappa:= elt_invariants(tmp_Phi, traces_tmp_Phi,kappa);
  until kappa.f = F;

  kappa.elt := Coerce(PR,kappa.elt);

  return kappa;

end;




#////////////////////////////////////////////////////////////////////////


#//-------------------------
round_four := function(Phi, with_certificates)
local PR,K,pi,Phihat,traces_Phi,traces_Phihat,P,S,P_v_star,E,phi_old_v_star,phi,F,kappa,nu,Pi,theta,
       s,Gamma,BETA,phis,poly,psi,poly_inv,gamma,beta,newPhi,r_L,rho_passes_hensel_test,delta,rho,         deltapsi,further_factorization;
#//-------------------------
#// the algorithm used is a hybrid between
#// pauli: factoring polynomials over local fields
#// and
#// ford-pauli-roblot: a fast algorithm for factoring polynomials over Qp
#//
#// step g) [extend the ground field] in the pauli algorithm is replaced by
#// by the combination of step 12 in algorithm 5.1 and step 2,3 in 5.2 in
#// ford-pauli-roblot
#//
#// if certificates is false the information needed for the certificates is not
#// computed



 #   vprint RoundFour,1:"R4: factoring Phi with deg(Phi) =",Degree(Phi);
 #  vprint RoundFour,3:"R4: Phi =",Phi;

 #   FACTOR := recformat
 #             <factor:RngUPolElt,
 #              multiplicity:IntegerRing(),
 #              F:IntegerRing(),              // inertia deg
 #              Gamma:RngUPolElt,          // unramified
 #              Gamma_denval:IntegerRing(),   // certificate
 #              E:IntegerRing(),              // ram index
 #              Pi:RngUPolElt,             // ramified
 #              Pi_denval:IntegerRing(),      // certificate
 #              IdealGen1,                 // two element representation
 #              IdealGen2                  // of the prime ideal
 #             >;

further_factorization := function(factors, with_certificates)
local L,i,certs,j,fact;
  L := [];
  for i in [1..Length(factors)] do
      fact:=factors[i][1];
      while IsList(fact) do
           fact:= fact[1];
      od;
      certs:=round_four(fact, with_certificates);
    for j in [1..Length(certs)] do
      Append_(L,[certs[j]]);
    od;
  od;
  return L;
end;

    PR := Parent(Phi);
    K := CoefficientRing(PR);
    pi := UniformizingElement(K);

  #  vprint RoundFour,4:"R4: over",K;

    if Degree(Phi) <= 1 then
      return [rec(factor:=Phi,multiplicity:= 1, F:=1,Gamma:=Coerce(PR,1),Gamma_denval:=0, E:=1,Pi:=Coerce(PR,pi),Pi_denval:=0)];
    fi;


    Phihat := Phi;
    traces_Phi := traces_from_poly(Phi);
    traces_Phihat := traces_Phi;
    #//////////////////////////// RAMIFIED DATA
    P := [];                  #// list of polynomials that allows
                              #// construction of polynomials with
                              #// v_star valuation with den up to E
    S := [];                  #// list of maximal exponents to be used
                              #// for the polynomials in P
    P_v_star := [];           #// list of v_star valuations of the elements
                              #// of P
    E := 1;                   #// highest known ramification index
    phi_old_v_star := 0;      #// old v* valuation of phi, used to check
                              #// that valuation of phi increases
    #////////////////////////////
    BETA := eltRec(PR.1,0);      #// we record the change of the reference
                              #// polynomial Phihat in BETA/pi^BETA_denval
    #////////////////////////////

#//////////////////////////////////////////////////////////////////////////





    phi := eltRec(PR.1,0);
    phi.chi := Phihat;
    phi:=elt_invariants(Phihat, traces_Phihat, phi);

    if phi.passes_newton_test=FALSE or phi.passes_hensel_test=FALSE then
        #vprint RoundFour,1:"R4: x fails Newton or Hensel test -- factoring";
        return further_factorization(newton_lift_factorization_chi2(Phi, phi), with_certificates);
    fi;

    #///////////////////// UNRAMIFIED DATA
    F     := phi.f;    #// highest known inertia degree
    kappa := phi;      #// generator of the unramified extension of degree F
    nu    := kappa.nu; #// generating polynomial of the unramified
    #///////////////////// extension of degree F

    phi   := eltRec(phi.nu,0);
    phi:=elt_invariants(Phihat, traces_Phihat,phi);

    #vprint RoundFour,1:"R4: initializing E to",E,"and F to",F;
    #vprint RoundFour,3:"R4: inertia element:",BETA`elt;

    if F = Degree(Phi) then
      #vprint RoundFour,1:"R4: F =",F,"= deg(Phi)  =>  Phi is irreducible";
      #//vprint RoundFour,1:"R4: computing 2 element certificate";
      return [rec(factor:=Phi,multiplicity:=1,F:=F,Gamma:=PR.1,Gamma_denval:=0,E:=1, Pi:=Coerce(PR,pi), Pi_denval:= 0)];
    fi;

    #//print "nu_phi",nu_phi;


    #// main loop
    while true do
       # vprint RoundFour,2:"R4| phi =",phi`elt;
       #vprint RoundFour,1:"R4| deg(phi) =", Degree(phi`elt),"and v*(phi) =", phi`v_star;
       # vprint RoundFour,2:"R4| E =",E,"and F =",F;
        #/* Step a ***************************************************************/
	#Print(phi.v_star," ",phi_old_v_star,"\n");
        if  phi.v_star <= phi_old_v_star then
          return Error ("Factorization: Insufficient precision in Round Four main loop");
        fi;
        if Degree(phi.elt) > E*F then
          return Error ("Factorization: Degree of phi too large in Round Four main loop");
        fi;

        if phi.passes_newton_test=FALSE then
	    #vprint RoundFour,1:"R4a phi fails Newton test -- factoring";
            theta := elt_evaluate_trans(Phi,traces_Phi,phi,BETA);
	    return further_factorization(newton_lift_factorization_chi2
                              (Phi, theta),with_certificates);
	    #vprint RoundFour,4:"R4a phi passes Newton test";
	fi;

	#/* Step b [increase E] **************************************************/
        if E mod phi.e <> 0 then
	    Append_(P, [phi.elt]);
	    Append_(P_v_star, [v_star(chi(Phihat, traces_Phihat, phi.elt))]);
	    s := Div(LCM(E, phi.e),E);
	    Append_(S, [s]);
            #vprint RoundFour,1:"R4b increasing E from",E,"to",s*E;
            #vprint RoundFour,4:"R4b list of polynomials giving ramification",P;
            #vprint RoundFour,4:"R4b with v* valuations",P_v_star;
            #vprint RoundFour,4:"R4b and maximal exponents",S;
	    E := s * E;

            if E*F = Degree(Phihat) then #// Phi is irreducible
              #vprint RoundFour,1:"R4b E * F =",E,"*",F,"= deg(Phi) =",Degree(Phihat),
               #                  "  =>   Phi is irreducible";
              if with_certificates then
                #vprint RoundFour,1:"R4b computing 2 element certificate";
                #//vprint RoundFour,5:"R4b computing ramified certificate";
                Pi := elt_with_v_star(Phihat, P, P_v_star, S, E, 1/E);
                Pi := elt_evaluate(Phi, Pi, BETA);
                Pi:=elt_invariants(Phi, traces_Phi, Pi);
                if Pi.v_star <> 1/E then
                  return Error ("Factorization: Insufficient precision in computation of ramified certificate");
                fi;
                #//vprint RoundFour,5:"R4b computing unramified certificate";
                if F = 1 then
                  Gamma := eltRec(Coerce(Parent(Phi),1),0);
                else
                  Gamma := BETA;
                fi;
                Gamma:=elt_invariants(Phi, traces_Phi, Gamma);
                #//print "F",F;
                if Gamma.f <> F then
                  return Error ("Factorization: Insufficient precision in computation of unramified certificate");
                fi;
                return [rec(factor:=Phi, multiplicity:=1,F:=F, Gamma:=Gamma.elt, Gamma_denval:=Gamma.denval, E:=E, Pi:=Pi.elt,Pi_denval:=Pi.denval)];
              else #// no certificates
                return [rec( factor:=Phi, multiplicity:= 1,F:=F,E:=E)];
              fi;

           else #// continue in main loop
                 #// increase valuation of phi, such that deg(phi)=E*F
              phis := eltRec((phi.elt)^s,0); phis.v_star := s*phi.v_star;
              phi := phis;
              #vprint RoundFour,2:"R4b phi =",phi`elt;
              #vprint RoundFour,1:"R4b deg(phi) =", Degree(phi`elt),
               #                  "and v*(phi) =", phi`v_star;
           fi;
	fi;

	#/* Step c **************************************************************/
        #vprint RoundFour,5:"R4c find psi with v*(psi) = v*(phi) =",phi`v_star;
        poly := elt_with_v_star(Phi, P, P_v_star, S, E, phi.v_star);
        psi := poly.elt * pi^-poly.denval;

        #/* Step d *************************************************************/
        #vprint RoundFour,5:"R4d gamma := phi/psi";
        poly_inv := elt_inverse(poly, Phihat);

        gamma := eltRec((phi.elt * poly_inv.elt) mod Phihat, poly_inv.denval-poly.denval);
        gamma.elt := gamma.elt mod pi^(gamma.denval+1);
        gamma:=elt_invariants(Phihat, traces_Phihat, gamma);

        if gamma.v_star <> 0 then
          #//print vstar_phi,old_vstar_phi;
          return Error ("Factorization: Insufficient precision in Round Four main loop");
        fi;

        #/* Step e *************************************************************/
	if gamma.passes_hensel_test=FALSE then
	    #vprint RoundFour,1:"R4e gamma fails Hensel test -- factoring";
            theta := elt_evaluate_trans(Phi,traces_Phi,gamma,BETA);
	    return further_factorization(hensel_lift_factorization_chi2
                                         (Phi, theta), with_certificates);
        #else
           #Print ("RoundFour,4:R4e gamma passes Hensel test");
	fi;

	#/* Step f ************************************************************/
        #//vprint RoundFour,5:"R4f new inertia by gamma ?";
        if F mod gamma.f <> 0 then #// increase F
          #// change reference polynomial Phihat, such that
          #//   F_x = LCM(F,F_gamma)
          #// start the main loop again with  Step a
          kappa := ff_composition_gen(Phihat, PR.1, F, gamma);
          if E * LCM(F, gamma.f) = Degree(Phihat) then #// Phi is irreducible
             F := LCM(F, gamma.f);
             #vprint RoundFour,1:"R4f E * lcm(F,F_gamma) =",E,"*",F,"= deg(Phi) =",Degree(Phihat),
              #               "  =>   Phi is irreducible";
             if with_certificates then
		#vprint RoundFour,1:"R4f computing 2 element certificate";
                Pi := elt_with_v_star(Phihat, P, P_v_star, S, E, 1/E);
                Pi := elt_evaluate(Phi, Pi, BETA);
                Pi:=elt_invariants(Phi, traces_Phi, Pi);
                if Pi.v_star <> 1/E then
                  return Error ("Factorization: Insufficient precision in computation of ramified certificate");
                fi;
                Gamma := elt_evaluate(Phi, kappa, BETA);
                Gamma:=elt_invariants(Phi, traces_Phi, Gamma);
                #//print "F",F;
                if Gamma.f <> F then
                  return Error ("Factorization: Insufficient precision in computation of unramified certificate");
                fi;

                return [rec(factor := Phi, multiplicity := 1,
                            F:=F, Gamma := Gamma.elt, Gamma_denval := Gamma.denval,
                            E:=E, Pi := Pi.elt, Pi_denval := Pi.denval)];
             else
               return [rec(factor := Phi, multiplicity := 1,F:=F,E:=E)];
             fi;
	  #/* Step g ************************************************************/
         else #// find a new reference polynomial Phihat
            F := LCM(F, gamma.f);
	    #vprint RoundFour,1:"R4 increasing F from",F,"to", F, "and resetting E to 1";
#//"kappa",kappa;
            beta := kappa;
            #//beta`elt := beta`elt+pi^(beta`denval+1)*PR.1;
            #vprint RoundFour,5:"R4g searching beta with F_beta = F =",F,
             #                  "and chi_beta irreducible of Degree deg(Phi) =",Degree(Phi);
            newPhi := elt_chi(Phihat, traces_Phihat, beta);
            #//while not is_squarefree(newPhi) do
            if Discriminant(Phihat) <> 0 then
              #// we can use the discriminant for testing whether newPhi is
              #// squarefree
              while Discriminant(newPhi)= 0 do
                #// assure that the new reference polynomial is squarefree
                #// then it generates the same extension as Phi
                #vprint RoundFour,5:"R4g beta := beta + pi*x";
                beta.elt := beta.elt+pi^(beta.denval+1)*PR.1;
                newPhi := elt_chi(Phihat, traces_Phihat, beta);
              od;
            else
              #// better also use the gcd for testing whether newPhi is
              #// squarefree as we have a precion problem with the discriminant
              #// already
              while is_squarefree(newPhi)=FALSE do
                #vprint RoundFour,5:"R4g beta := beta + pi*x";
                beta.elt := beta.elt+pi^(beta.denval+1)*PR.1;
                newPhi := elt_chi(Phihat, traces_Phihat, beta);
              od;
            fi;

            Phihat := newPhi;
            traces_Phihat := traces_from_poly(Phihat);

            #// reinitialize
            E := 1;
            P := [];
            S := [];
            P_v_star := [];
            phi_old_v_star := 0;

            BETA := elt_evaluate(Phi, beta, BETA);
            #vprint RoundFour,4:"R4g new inertia element beta =",BETA`elt;
            #vprint RoundFour,4:"R4g with denominator pi ^",BETA`denval;

            phi := eltRec(PR.1,0);
	    phi.chi := Phihat;
            phi:=elt_invariants(Phihat,traces_Phihat,phi);
            if (phi.passes_newton_test and phi.passes_hensel_test) =FALSE then
              #vprint RoundFour,1:"R4g Newton test or Hensel test failed -- factoring";
              theta := elt_evaluate_trans(Phi,traces_Phi,phi,BETA);
              return further_factorization(newton_lift_factorization_chi2
                       (Phi, theta), with_certificates);
            fi;

            nu := phi.nu;
            phi:= eltRec(phi.nu,0);
	    phi:=elt_invariants(Phihat,traces_Phihat,phi);
	  fi;
	else
          #// no increase of F
          #// continue in main loop
        #/* Step h *************************************************************/
         # vprint RoundFour,5:
          #"R4h find integral delta with chi(delta)=chi(gamma) mod (pi) and deg(delta)<F";
          r_L:= res_field_rep(Phihat,traces_Phihat,nu,gamma);
	  rho_passes_hensel_test:=r_L[1];
	  delta:= r_L[2];
	  rho:= r_L[3];
          if rho_passes_hensel_test=FALSE then
	    #vprint RoundFour,1:"R4h rho fails Hensel test -- factoring";
            theta := elt_evaluate_trans(Phi,traces_Phi,rho,BETA);
	    return further_factorization(hensel_lift_factorization_chi2
                   (Phi, theta),with_certificates);
          fi;
          #vprint RoundFour,3:"R4h delta =",delta;
        #/* Step i *************************************************************/
         # vprint RoundFour,5:"R4i phi := phi - delta * psi";
	  deltapsi := delta * psi;
          #//print "deg", Degree(deltapsi);
          phi_old_v_star := phi.v_star;
          phi := eltRec(phi.elt-deltapsi,0);
	  phi:=elt_invariants(Phihat,traces_Phihat,phi);
        fi; #// increase F
    od;
end;



#////////////////////////////////////////////////////////
#// auxilary functions for factorization.

#// return true, g, dival, uni, exp such that
#// g := uni*pi^multval*f(pi^exp*x) has integral coefficients and is monic
#// g is over the quotient ring !
#// if because of insuffient precision such elements cannot be computed, return
#// false, 0, 0, required precision, exp

#  // if f is monic and integral do nothing
#  //if IsMonic(f) and Minimum([Valuation(a) : a in Eltseq(f)]) ge 0 then
#  //
#  //return

make_monic_and_integral := function(f)
local PK,K,overfield,R,PR,elseq,vallist,minval,divi,g,multval,PR,pi,lead,vallead,exp,expmult,gList,prec,QR,PQR,uni,SList;
  PK := Parent(f);
  K  := CoefficientRing(PK);
  overfield := IsField(K);
  if overfield then
    R := IntegerRing(K);
    PR := PolynomialAlgebra(R);
    elseq:= ElementToSequence(f);
    vallist:= [1..Length(elseq)];
    Apply_(vallist,x->Valuation(elseq[x]));
    minval := Minimum(vallist);
    divi := UniformizingElement(K)^(minval);
    g := Coerce(PR,(f/divi));
    multval := -minval;
  else
    R := K;
    PR := PK;
    elseq:= ElementToSequence(f);
    vallist:=[1..Length(elseq)];
    Apply_(vallist,x->Valuation(elseq[x]));
    minval := Minimum(vallist);
    divi := UniformizingElement(K)^(minval);
    g := Coerce(PR,(Div(f,divi)));
    multval := -minval;
  fi;
  #//print g;
  PR := PolynomialAlgebra(R);
  pi := UniformizingElement(R);
  lead := LeadingCoefficient(g);
  vallead := Valuation(lead);
  elseq := ElementToSequence(g);
  exp := -1;
  if Degree(g) > 0 then
    repeat
      exp :=exp+1;
      expmult := Degree(g)*exp-vallead;
      #//print "expmult",expmult;
      vallist:= [2..Length(elseq)];
      Apply_(vallist,i->Valuation(elseq[i])-exp*(i-1)+expmult);
      minval:= Minimum(vallist);
    until (expmult >= 0 and minval >= 0);
    if vallead+expmult >= Precision(R) or exp*Degree(g) >= Precision(R) then
      return [false, 0, 0, vallead+expmult+exp*Degree(g), exp];
    fi;
  else
    exp := 0;
    expmult := - vallead;
  fi;
  g := g * pi^expmult;
  #//print "exp",exp;
  #//print g;
  SList:= ElementToSequence(g);
  gList:=[1..(Degree(g)+1)];
  Apply_(gList,i->Div(SList[i],pi^(exp*(i-1))));
  #//print g;
  multval := multval+expmult;
  if not Is (Type(R),res^pad) and not Is (Type(R),res^ser) then
     if Precision(R)=INFTY then prec:= _realPrec;
     else prec:= Precision(R);
     fi;
     if Is(Type(R),ord^ser) or Is(Type(R),ord^pow) then
        QR := Quotient(PowerSeriesRing(CoefficientRing(R)),prec);
     else
     QR := Quotient(R,pi^prec);
     fi;
  else
    QR := R;
  fi;
  PQR:= PolynomialAlgebra(QR);
  g := Coerce(PQR,g);
  uni := LeadingCoefficient(g);
  #//print uni;
  g := Div(g,uni);
  if IsMonic(g)=FALSE then return Error("g not monic"); fi;
  return [true, g, uni, multval, exp];
end;

#// make the reverse changes made to f by make_monic_and_integral to h
#// actually only the exp part, the rest goes into the scalar returned
#// with the polynomial factorization.
undo_monic_and_integral := function(f,exp,h)
local g,PK,K,pi,g;
  g := h;#//*uni;

  PK := Parent(f);
  K := CoefficientRing(PK);
  pi := UniformizingElement(K);

  g := Coerce(PK,g);
  g := Evaluate(g,pi^exp*PK.1);

  return g;
end;



#/////////////////////////////////////////////////////////////////////

#{For a polynomial f over a local ring or field return a precision at which
# the factorization of f as given by Factorization(f) will be Hensel liftable
# to the correct factorization.}


SuggestedPrecision:= function(f)
local max,seq,vallist,tmp_R,tmpi,mo,ok,h,uni,multval,exp,mmai_prec,PK,PR,P,valList,newloss,loss,gcd,
      old_P,dh,old_gcd,pgcd,done,good,hh,val;
  #// the trivial cases
  if Degree(f) = 0 then
    max := Maximum([2,Valuation(TrailingCoefficient(f))+1]);
    return max;
  elif Degree(f) = -1 then
    return 1;
  fi;

  #// monic and integral
  seq:= ElementToSequence(f);
  vallist:=[1..Length(seq)];
  Apply_(vallist,x->Valuation(seq[x]));
  if ((IsUnit(LeadingCoefficient(f)) and
       not IsField(CoefficientRing(Parent(f)))) or
      IsMonic(f)) and
      Minimum(vallist) >= 0 then
    #// make sure it is also squarefree
    if Valuation(Discriminant(f)) <
      Precision(CoefficientRing(Parent(f)))
      and Precision(CoefficientRing(Parent(f))) < INFTY then
      max := Maximum([2*RamificationDegree(CoefficientRing(Parent(f))),
                      2*Valuation(Discriminant(f))]);
      return max;
    fi;
  fi;
  tmp_R := CoefficientRing(Parent(f));
  tmpi:=UniformizingElement(tmp_R);
  if Is(Type(tmp_R),ord^ser) or Is(Type(tmp_R),res^ser) then
     f := Coerce(PolynomialAlgebra( Quotient(tmp_R,(2*Precision(tmp_R)))),f);
  else
     f := Coerce(PolynomialAlgebra( Quotient(tmp_R,tmpi^(2*Precision(tmpi)))) ,f);
  fi;
  mo:=make_monic_and_integral(f);
  ok:= mo[1];
  h:= mo[2];
  uni:=mo[3];
  multval:= mo[4];
  exp:= mo[5];
  mmai_prec := Valuation(LeadingCoefficient(f))+multval+exp*Degree(f);
  if not ok then
    #//print "increasing prec";
    mmai_prec := multval+exp*Degree(f);
    PK := PolynomialAlgebra(ChangePrecision(CoefficientRing(Parent(f)),
                                         mmai_prec));
    mo :=make_monic_and_integral(Coerce(PK,f));
    good:= mo[1];
    h:= mo[2];
    uni:= mo[3];
    multval:= mo[4];
    exp:= mo[5];
    if not good then
      return Error ("I cannot suggest a precision.");
    fi;
  fi;

  PR := Parent(h);
  P := CoefficientRing(Parent(h));

  #//print h;

   # // correct the precision for the derivative
    seq := ElementToSequence(h);
    valList:=[2..Length(seq)];
    Apply_(valList,x->Valuation(seq[x])+Valuation(Coerce(P,(x-1))));
    max := Maximum(valList);
    #//print "precision increase derivative",max-Precision(P);
    newloss := Maximum([max-Precision(P),0]);
#//RamificationDegree(P) eq 1 and InertiaDegree(P) eq 1 then
    repeat #// loop for the valuation of the discriminant
       loss := 0;
       gcd := 0;
       repeat
         old_P := P;
	 if not Precision(P) =INFTY then
         P := ChangePrecision(P,Precision(P)+(newloss));
	 fi;
         #//print "precision",Precision(P);
         PK := PolynomialAlgebra(ChangePrecision(CoefficientRing(Parent(f)),
                                         Precision(P)));
	 mo:=make_monic_and_integral(Coerce(PK,f));
         good:= mo[1];
         h:= mo[2];
         uni:=mo[3];
         multval:= mo[4];
         exp:= mo[5];
         #//print "new h",h;
         if not good then
            return Error ("I cannot suggest a precision.");
         fi;
         h := Coerce(PolynomialAlgebra(P),h);
         dh := Derivative(h);
         old_gcd := gcd;
         #//print "h",h;
         #//print "dh",dh;
         #vprint  RoundFour,4: "R4: gcd";
	 pgcd:=PGcd(h,dh);
         gcd:= pgcd[1]; newloss := pgcd[2];
         #//print "newloss",newloss;
         loss := loss + newloss;
         #//print "gcd",gcd, newloss;
         #//print "h mod gcd",h mod gcd;
         #//print "h div gcd",h div gcd;
       until (old_gcd = Coerce(PolynomialAlgebra(old_P),gcd) and h mod gcd = 0);
       hh := Div(h,gcd);
       #//print hh;
       if RamificationDegree(P) = 1 and InertiaDegree(P) = 1 then
         val := Valuation(Discriminant(Coerce(PolynomialAlgebra(RationalField()),hh)),
                          Prime(P));
         done := val.base <> INFTY;
         #//print val;
       else
         val := Valuation(Discriminant(Coerce(PolynomialAlgebra(P),hh)));
         done := val < Precision(P);
       fi;
       #//print "val",2*val,loss,mmai_prec;
    until done and newloss < Precision(P);
    max := Maximum([2*val+loss,mmai_prec+loss,2,Precision(P)]);
    return max;
end;

#{Factorization of a polynomial f over a local ring or a local field into monic irreducible
# factors together with the appropriate scalar and certificates for the irreducibility of the
# factors if requested }

LocalFactorization:= function(f,Certificates,IsSquarefree,Ideals,Extensions)
local with_certificates,mo,ok,h,uni,multval,exp,PQR,QR,sqrfr,is_sq_free,sq_free_part,loss,r4,res1,g,        i,hh,rem,zero,done,sequ,fa,valList,a,j,PR,R,pi,min,sca,K,PK,Kpi,cert,poly,copoly,gcd,G,C,D,Phi,        traces_Phi,Rho,U,PU,Pi,tmp_U,g2,tmp_PU,fact,T,eltseq,vallist,res2,qr,prec,un,tmpPR,tmpf;
    with_certificates := Certificates or Ideals or Extensions;

    if Precision(CoefficientRing(Parent(f))) = INFTY and
       Extensions then
      return Error ("Factorization: 'Extensions' can only be used over finite precision rings/fields");
    fi;

    #//print "Factorization",f;
    mo:=make_monic_and_integral(f);
    ok:= mo[1];
    h:=mo[2];
    uni:=mo[3];
    multval:=mo[4];
    exp := mo[5];
    if not ok then
      return Error ("Factorization: Insufficient precision");
    fi;
    #//print h;
    PQR := Parent(h);
    QR := CoefficientRing(PQR);

    sqrfr:= is_squarefree(h);
    is_sq_free:= sqrfr[1];
    sq_free_part:=sqrfr[2];
    loss:=sqrfr[3];
    if is_sq_free then
     
      r4 := round_four(h,with_certificates);
      res1:=[1..Length(r4)];
      Apply_(res1,x->Tuple([undo_monic_and_integral(f, exp, r4[x].factor),1]));
    else
      if loss>0 and not Is(Type(f),elt-alg^pol/ord^ser) or Is(Type(f),elt-alg^pol/res^ser) then
         if Precision(QR) <>INFTY then
            tmpPR:= _ChangePrecisionCoeff(PQR,Precision(QR)+loss);
         else 
            tmpPR:= _ChangePrecisionCoeff(PQR,_realPrec+loss);
         fi;
         tmpf:= Coerce(tmpPR,f);
         h:= make_monic_and_integral(tmpf)[2];
         sqrfr:= is_squarefree(h);
         sq_free_part:=sqrfr[2];
      fi;
      g := sq_free_part;
      fa:= 1;
#//"h",h;
#//"g",g;
      while IsMonic(g)=false do
         # Print("no mo\n");
         #mo:= make_monic_and_integral(g)[2];
         fa:= fa*LeadingCoefficient(g);
         mo := g*(IsInvertible(LeadingCoefficient(g)).ext1);
         g:= is_squarefree(mo)[2];
      od;
      r4 := round_four(g,with_certificates);
      if loss>0 and not Is(Type(f),elt-alg^pol/ord^ser) or Is(Type(f),elt-alg^pol/res^ser) then
         for i in [1..Length(r4)[1]] do
             if Precision(QR)<>INFTY then
                r4[1][i].factor:= Coerce(PQR,r4[1][i].factor);
             else tmpPR:= _ChangePrecisionCoeff(PQR,_realPrec);
                r4[1][i].factor:= Coerce(tmpPR,r4[1][i].factor);
                r4[1][i].factor:= Coerce(PQR,r4[1][i].factor);
             fi;
         od;
      fi;
      res1:=[1..Length(r4)];
      Apply_(res1,x->Tuple([r4[x].factor,0]));
      #//hh := h;
      for i in[1..Length(res1)] do
        hh := h;
        repeat
          #//print "hh pre",hh;
          hh := Div(hh,res1[i][1]);
          rem := hh mod res1[i][1];
#Print(" = hh mod res ",rem,"\n");
          res1[i][2] := res1[i][2]+1;
          #//print "hh post",hh;
          zero := (rem = Coerce(Parent(rem),0));
          done := false;
          if not zero then
	    sequ:= ElementToSequence(rem);
	    valList:=[1..Length(sequ)];
	    Apply_(valList,x->Valuation(sequ[x]));
            done := Minimum(valList) < Precision(QR)-loss;
          fi;
        until done;
      od;
      #//print res;
#//print "exp sum", &+[a[2]*Degree(a[1]) : a in res],"=?=",Degree(f);
 #     //print res;
       a:= res1[1][2]*Degree(res1[1][1]);
       for j in [2..Length(res1)] do
           a:= a+res1[j][2]*Degree(res1[j][1]);
       od;
      if a<> Degree(f) then
        Print("a=",a);
#//"res",res;
#//"f",f;
 return Error ("Factorization: Insufficient precision in factorization of a non-squarefree polynomial.");
      fi;
      res2:= [1..Length(res1)];
      Apply_(res2,x->Tuple([undo_monic_and_integral(f, exp, res1[x][1]), res1[x][2]]));
      res1:= res2;
    #//print res;
    fi; #// squarefee by gcd
#//  end if;  // is_squarefree_easy

    PR := Parent(f);
    R := CoefficientRing(PR);
    pi := UniformizingElement(R);
    if not IsField(R) then  #// distribute the scalar factor 1/pi^multval
      for i in [1..Length(res1)] do #// over the factors
        eltseq:= ElementToSequence(res1[i][1]);
	vallist:=[1..Length(eltseq)];
	Apply_(vallist,x->Valuation(eltseq[x]));
        min :=  Minimum(vallist);
        res1[i][1]:=Div( res1[i][1],pi^min);
        multval := multval -res1[i][2]*min;
      od;
      if multval > 0 then
        return Error ("Factorization: Something went wrong when trying to factor a non-monic Polynomial.");
      fi;
    fi;
    sca := Coerce(R,uni) * pi^(-multval);

   # //==========================================
    if with_certificates then
   #   CERT := recformat<F:IntegerRing(),   // inertia deg
   #                     Rho:RngUPolElt, // unramified certificate
   #                     E:IntegerRing(),   // ram index
   #                     Pi:RngUPolElt,  // ramified certificate
   #                     IdealGen1,      // 2 element representation of
   #                     IdealGen2:RngUPolElt,// the corresponding global ideal
   #                     Extension>;     // extension of R given by the
                                    #    // factor
      if not IsField(R) then
        #if Is(Type(R),res^pad) then
	   #K:= pAdicField(Coerce(Z,UniformizingElement(R)),Precision(R));
	#else
           K := FieldOfFractions(R);
	#fi;
        PK := PolynomialAlgebra(K);
      else
        K := R;
        PK := PR;
      fi;


      Kpi := UniformizingElement(K);
      cert:=[1..Length(r4)];
      Apply_(cert,x->rec(F:= r4[x].F, Rho:=Coerce(PK,((r4[x].Gamma))/Kpi^r4[x].Gamma_denval),
                         E:=r4[x].E,  Pi:=Coerce(PK,Coerce(PK,r4[x].Pi)/Kpi^r4[x].Pi_denval)));

      #//-----------------------------------------
      if Ideals then
        for i in [1..Length(res1)] do
          cert[i].IdealGen1 := pi;
          poly := Coerce(PK,res1[i][1]);
          copoly := 1;
          for j in [1..Length(res1)] do
            if j <> i then
              copoly := copoly *Coerce(PK,res1[j][1]);
            fi;
          od;
	  gcd:=XGCD(poly,copoly);
          G := gcd.base;
	  C:= gcd.ext1;
	  D:= gcd.ext2;
          cert[i].IdealGen2 := cert[i].Pi*D*copoly+poly*C;
        od;
      fi;
      #//-----------------------------------------
      if Extensions then
        #vprint RoundFour,1:" -- computing local extensions";
        for i in [1..Length(res1)] do
            Phi := r4[i].factor;
            #//print "Phi",Phi;
            traces_Phi := traces_from_poly(Phi);
            if cert[i].F > 1 then
              Rho := eltRec(r4[i].Gamma, r4[i].Gamma_denval);
              Rho:=elt_invariants(Phi, traces_Phi, Rho);
             # //print "Gamma",r4[i]`Gamma, r4[i]`Gamma_denval;
             # //print "chi_Rho",chi_Rho;
              if not Rho.passes_hensel_test then
               return Error ("Factorization: something went wrong with the unramified certificate");
              fi;
              if Rho.f <> cert[i].F then return Error("Error in Extensions");fi;
              U := UnramifiedExtension(R,Rho.nu);
            else
              U := R;
            fi;
            if cert[i].E > 1 then
              Pi := eltRec(r4[i].Pi, r4[i].Pi_denval);
              Pi:=elt_invariants(Phi, traces_Phi, Pi);
              PU := PolynomialAlgebra(U);
              #// unfortunately the precison used to factor f must noty be
              #// sufficient for factoring chi_Pi
              tmp_U := U;
              #//tmp_U := ChangePrecision(U,
              #//    Maximum([2*Valuation(Discriminant(chi_Pi)),Precision(R)]));
              tmp_PU := PolynomialAlgebra(tmp_U);
              #//vtime RoundFour,4:
              fact := LocalFactorization(Coerce(tmp_PU,(Coerce(PU,(Coerce(PR,Pi.chi))))),false,false,false,false);
              #//"fact4",fact;
	      if Degree(fact[1][1][1]) <> Div(Degree(Phi),cert[i].F) then
	         return Error("Degree");
	      fi;
              if (2*Valuation(Discriminant(fact[1][1][1]))+cert[i].E)/cert[i].E
                 > Precision(U) then
                  return Error ("Factorization: The Precision was not high enough",
                       "to assure that the ramified extension is unique");
              fi;
              T := TotallyRamifiedExtension(U,Coerce(PU,fact[1][1][1]));
            else
              T := U;
            fi;
            cert[i].Extension := T;
        od;
      fi;
      return [res1,sca, cert];
    #//===========================================
    else
      return [res1,sca];
    fi;
end;



#{True iff f is irreducible}
IsIrreducible:= function(f)
local F;
    F := Factorization(f);
    return Length(F) = 1;
end;


#{Factorization of a polynomial f over a local ring or a local field into monic irreducible
# factors together with the appropriate scalar and certifcates for the irreducibility of the
# factors if requested }

#Factorization:= function(f,Certificates,IsSquarefree,Ideals,Extensions)
#    return LocalFactorization(
#	f ,Certificates, IsSquarefree ,Ideals,Extensions);
#end;

InstallMethod(
rec(
name := "Factorization",
kind := "FUNCTION",
sin := [[elt-alg^pol/any^loc,"Phi"]],
sou := [[list]],
opt := [[elt-alg^boo,"Certificates",
         "If TRUE two element certificates for the irreducibilty of factors are returned.", rec(Default:=FALSE)],
        [elt-alg^boo,"IsSquarefree",
         "If TRUE the polynomial is assumed to be squarefree",rec(Default:=FALSE)],
        [elt-alg^boo,"Ideals","",rec(Default:=FALSE)],
        [elt-alg^boo,"Extensions",
         "If TRUE the extensions generated by the irreducible factors are returned",
	 rec(Default:=FALSE)]],
short := "The factorization of a polynomial over a p-adic field or a "+
" p-adic ring.  The algorithm used is a combinantion of the algorithms by "+
"Ford-Pauli-Roblot and by Pauli."
),
function(arg)
  local opt;
#Print("arg",arg,"\n");
  opt := Optarg(arg,
                ["Certificates", "IsSquarefree", "Ideals", "Extensions"],
                [FALSE,FALSE,FALSE,FALSE]);
   
#Print("opt",opt,"\n");
#Print("call",arg[1],opt.Certificates, 
#                              opt.IsSquarefree, 
#                              opt.Ideals,
#                              opt.Extensions,"\n");
  return LocalFactorization(arg[1],opt.("Certificates"), 
                              opt.("IsSquarefree"), 
                              opt.("Ideals"),
                              opt.("Extensions"));
#  return LocalFactorization(arg[1],TRUE,FALSE,FALSE,FALSE);
end
);

#q3 := pAdicRing(3,30);
#q3x := PolynomialRing(q3);
#f := (q3x.1^3+3)*(q3x.1^3+6);
#LocalFactorization(f,FALSE,FALSE,FALSE,FALSE);
#ReadLib("locFact");
#Factorization(f);
#q3 := pAdicRing(3,30);
#q3x := PolynomialRing(q3);

