// ****************************** MODULAR ******************************
zx<x>:=PolynomialRing(Integers());
print "zx is", zx;

// FactMod(f,p) returns the factorization of the integer polynomial f over Fp.

FactMod := function(f,p);
  return Factorization(PolynomialRing(Integers(p))!f), 
  PolynomialRing(Integers(p))!f;
end function;

// facdeg(f,p) returns the partition of (Degree(f)) in the degrees of the
// irreducible factors of f modulo p, in the following way:
// facdeg(f,p) = [i_1,i_2,...,i_n], means that there are i_j irreducible factors
// of degree j.

facdeg := function(f,p);
  if IsDivisibleBy(LeadingCoefficient(f),p) then
    return [Degree(f)];
  end if;
  fact := Factorization(PolynomialRing(Integers(p))!f);
  facdeg := zeroseq(Degree(fact[#fact][1]));
  for i := 1 to #fact do
    facdeg[Degree(fact[i][1])] +:= fact[i][2];
  end for;
  return facdeg;
end function;

// fills(fd1,fd2) return whether fd1 is a subpartition of fd2, or not.

fills := function(fd1,fd2);
  if #fd2 le 1 then return true; end if;
  gb := #fd2;
  fd2[#fd2] -:= 1;
  if fd2[#fd2] eq 0 then fd2 := removezeros(fd2); end if;
  tries := [k:k in [gb..#fd1] | fd1[k] ne 0];
  if #tries eq 0 then return false; end if;
  for try in tries do
    newfd1 := minusbofk(fd1,gb,try);
    if $$(newfd1,fd2) then return true; end if;
  end for;
  return false;
end function;

// IsParentOf(fd1,fd2) returns whether fd1 is a parent of fd2.
// First there are some simple rules that might determine this,
// if this does not work "fills" is used.

IsParentOf := function(fd1,fd2);
  if #fd1 lt #fd2 then return false; end if;
  
  // Next we remove equal parts.
  for i := 1 to #fd2 do
    m := Min(fd1[i],fd2[i]);
    if m ne 0 then
      fd1[i] -:= m; fd2[i] -:= m;
    end if;
  end for;
  
  fd1 := removezeros(fd1); fd2 := removezeros(fd2);
    
  fd1sum := 0; fd2sum := 0;
  for c := 1 to #fd2 do
    fd1sum +:= c*fd1[c];
    fd2sum +:= c*fd2[c];
    if fd1sum gt fd2sum then return false; end if;
  end for;
  
  return fills(fd1,fd2);
end function;

// DirectParents(fds) returns the partitions that are direct parent of all 
// partitions in fds, this means a parent by taking adding just two elements

DirectParents := function(fds);
  // all fds have the same fdsum (we can assume that)
  dparents := {};
  for fd in fds do
    for a in {k:k in [1..#fd] | fd[k] ne 0} do
      if fd[a] eq 1 then initb := a+1;
        else initb := a;
      end if;
      for b in {l:l in [initb..#fd] | fd[l] ne 0} do
          dparents join:= {plusab(fd,a,b)};
      end for;
    end for;
  end for;
  return dparents;
end function;

// Parents(fd) gives all the parents of the partition fd, using "DirectParents"

Parents := function(fd);
  parents := {fd};
  dparents := {fd};
  
  while #dparents ne 0 do
    dparents := DirectParents(dparents);
    parents join:= dparents;
  end while;
  
  return parents;
end function;
  

// modcomb(f,B) tries to find a combination of primes =< B that give a
// `modulo p-certificate'. f has to be of degree >= 2.
      
modcomb := function(f,B);
  p := 2;
  hp := p;
  fd := facdeg(f,p); 
  mhf := &+fd;

  if #fd ne 1 then
    ps := [2]; facdegs := [fd]; possible := Parents(fd); beginpfound := true;
  else ps := []; facdegs := []; beginpfound := false;
  end if;
  
  while not beginpfound and p lt B do
    p := NextPrime(p);
    fd := facdeg(f,p);
    if #fd ne 1 then
      ps := [p]; facdegs := [fd]; possible := Parents(fd); beginpfound := true;
      hf := &+fd;
      if hf le mhf then
        mhf := hf;
        hp := p;
      end if;
    end if;
  end while;
  
  if p gt B then ps := [PreviousPrime(p)]; facdegs:= [facdeg(f,ps[1])];
                 possible:=[];
  end if;
  
  while #possible gt 1 and p lt B do
    p := NextPrime(p);
    fd := facdeg(f,p);
    if #fd eq Degree(f) then return {fd},[p],[fd],p; end if;
    if #fd ne 1 then 
      newpos := {};
      for poss in possible do
        if IsParentOf(poss,fd) then newpos := newpos join {poss}; end if;
      end for;
      // adapt Hensel prime hp if improvement occurs:
      hf := &+fd;
      if hf le mhf then
        mhf := hf;
        hp := p;
      end if;
      if #newpos lt #possible then
        Append(~ps,p); Append(~facdegs,fd); possible := newpos;
      end if;
    end if;
  end while;
  return possible,ps,facdegs,hp;
end function;

// facttookfact(fact) is an auxiliary algorithm to rewrite a factorization

facttookfact := function(fact);
  okfact := [];
  for i:=1 to #fact do
    for j:=1 to fact[i,2] do
      Append(~okfact,fact[i,1]);
    end for;
  end for;
  return okfact;
end function;

// fabs(f) returns |f| as needed for the calculation of the Mignotte bound

fabs := function(f);
  return Sqrt(&+[i^2 : i in Coefficients(f)]);
end function;

// MignotteBound(f) returns a bound on the coefficients of a non-trivial factor
// of f, based on the idea that the factor has degree =< Degree(f) / 2, 
// and on a bound by Mignotte.

MignotteBound := function(f);
  n := Degree(f);
  d := Floor(n/2);
  MB := Binomial(d-1,Floor(d/2))*fabs(f) + Binomial(d-1,Floor(d/2)-1);
  return MB;
end function;

// modular(f) tries to find a `modulo p-certificate', if such a certificate can
// not be found a `Hensel-certificate' is computed.
// input: integer polynomial f
// output: <i, S> where the possibilities for i are
//   i = 0, meaning that a true factor is found; the factor is S
//   i = 1, meaning that a modular certificate is generated (and returned as S;
//          as a 2-tuple of primes and the modular degrees of f mod primes)
//   i = 2, meaning that a Hensel certificate was found; S is the prime&power

modular := function(f);
  d := Degree(f);
  poss,ps,facdegs,hp := modcomb(f,30*d);
  if #poss eq 1 then return <1,<ps,facdegs>>; end if;
  
//  If we haven't found a factor after 30*d tries we know by testing
//  that it is sensible to use the Hensel-method.
  
  MB := MignotteBound(f);
/*
  min, minposi := Minimum([fdsum(fd) : fd in facdegs]);
  p := ps[minposi];
*/
  p := hp;
  min := &+facdeg(f,p);
  pow := Ceiling(Log(p,2*MB));
  P<w> := PolynomialRing(Integers(p));
  Q<z> := PolynomialRing(Integers(p^pow));
  pfact := Factorization(P!f);
  okpfact := facttookfact(pfact);
  
  henselfact := HenselLift(f,okpfact,Q);
  
  okparts := [part:part in Partitions(min) | #part eq 2];
  for okpart in okparts do
    tries := Subsets({1..min},okpart[2]);
    for try in tries do 
      g := 1;
      for i in try do
        g *:= zx!henselfact[i];
      end for;
      
      coefg := Coefficients(g);
      for i:=1 to Degree(g) do
        if coefg[i] gt (p^pow)/2 then g -:= (p^pow)*x^(i-1); end if;
      end for;
      
      if IsDivisibleBy(f,g) then return <0,g>; end if;
    end for;
  end for;
  return <2,p,pow>;
end function;