#include <NTL/ZZXFactoring.h>#include <NTL/lzz_pXFactoring.h>#include <NTL/vec_vec_long.h>#include <NTL/vec_vec_ulong.h>#include <NTL/vec_double.h>#include <NTL/new.h>NTL_START_IMPLstruct LocalInfoT {   long n;   long NumPrimes;   long NumFactors;   vec_long p;   vec_vec_long pattern;   ZZ PossibleDegrees;   PrimeSeq s;};staticvoid mul(ZZ_pX& x, vec_ZZ_pX& a)// this performs multiplications in close-to-optimal order,// and kills a in the process{   long n = a.length();   // first, deal with some trivial cases   if (n == 0) {      set(x);      a.kill();      return;   }   else if (n == 1) {      x = a[0];      a.kill();      return;   }   long i, j;   // assume n > 1 and all a[i]'s are nonzero   // sort into non-increasing degrees   for (i = 1; i <= n - 1; i++)      for (j = 0; j <= n - i - 1; j++)         if (deg(a[j]) < deg(a[j+1]))            swap(a[j], a[j+1]);   ZZ_pX g;   while (n > 1) {      // replace smallest two poly's by their product      mul(g, a[n-2], a[n-1]);      a[n-2].kill();      a[n-1].kill();      swap(g, a[n-2]);      n--;      // re-establish order      i = n-1;      while (i > 0 && deg(a[i-1]) < deg(a[i])) {         swap(a[i-1], a[i]);         i--;      }   }   x = a[0];   a[0].kill();   a.SetLength(0);}void mul(ZZX& x, const vec_pair_ZZX_long& a){   long l = a.length();   ZZX res;   long i, j;   set(res);   for (i = 0; i < l; i++)      for (j = 0; j < a[i].b; j++)         mul(res, res, a[i].a);   x = res;}void SquareFreeDecomp(vec_pair_ZZX_long& u, const ZZX& ff)// input is primitive {   ZZX f = ff;   ZZX d, v, w, s, t1;   long i;   u.SetLength(0);   if (deg(f) <= 0)      return;   diff(t1, f);   GCD(d, f, t1);   if (deg(d) == 0) {      append(u, cons(f, 1));      return;   }   divide(v, f, d);    divide(w, t1, d);   i = 0;   for (;;) {      i = i + 1;      diff(t1, v);      sub(s, w, t1);      if (IsZero(s)) {         if (deg(v) != 0) append(u, cons(v, i));         return;      }      GCD(d, v, s);      divide(v, v, d);      divide(w, s, d);      if (deg(d) != 0) append(u, cons(d, i));   }}staticvoid HenselLift(ZZX& Gout, ZZX& Hout, ZZX& Aout, ZZX& Bout,                const ZZX& f, const ZZX& g, const ZZX& h,                const ZZX& a, const ZZX& b, const ZZ& p) {   ZZX c, g1, h1, G, H, A, B;   mul(c, g, h);   sub(c, f, c);   if (!divide(c, c, p))      Error("inexact division");   ZZ_pX cc, gg, hh, aa, bb, tt, gg1, hh1;   conv(cc, c);   conv(gg, g);   conv(hh, h);   conv(aa, a);   conv(bb, b);   ZZ_pXModulus GG;   ZZ_pXModulus HH;   build(GG, gg);   build(HH, hh);   ZZ_pXMultiplier AA;   ZZ_pXMultiplier BB;   build(AA, aa, HH);   build(BB, bb, GG);   rem(gg1, cc, GG);   MulMod(gg1, gg1, BB, GG);      rem(hh1, cc, HH);   MulMod(hh1, hh1, AA, HH);   conv(g1, gg1);   mul(g1, g1, p);   add(G, g, g1);   conv(h1, hh1);   mul(h1, h1, p);   add(H, h, h1);   /* lift inverses */   ZZX t1, t2, r;   mul(t1, a, G);   mul(t2, b, H);   add(t1, t1, t2);   add(t1, t1, -1);   negate(t1, t1);   if (!divide(r, t1, p))      Error("inexact division");   ZZ_pX rr, aa1, bb1;   conv(rr, r);      rem(aa1, rr, HH);   MulMod(aa1, aa1, AA, HH);   rem(bb1, rr, GG);   MulMod(bb1, bb1, BB, GG);   ZZX a1, b1;   conv(a1, aa1);   mul(a1, a1, p);   add(A, a, a1);   conv(b1, bb1);   mul(b1, b1, p);   add(B, b, b1);   Gout = G;   Hout = H;   Aout = A;   Bout = B;}staticvoid HenselLift1(ZZX& Gout, ZZX& Hout,                 const ZZX& f, const ZZX& g, const ZZX& h,                const ZZX& a, const ZZX& b, const ZZ& p) {   ZZX c, g1, h1, G, H;   mul(c, g, h);   sub(c, f, c);   if (!divide(c, c, p))      Error("inexact division");   ZZ_pX cc, gg, hh, aa, bb, tt, gg1, hh1;   conv(cc, c);   conv(gg, g);   conv(hh, h);   conv(aa, a);   conv(bb, b);   ZZ_pXModulus GG;   ZZ_pXModulus HH;   build(GG, gg);   build(HH, hh);   rem(gg1, cc, GG);   MulMod(gg1, gg1, bb, GG);      rem(hh1, cc, HH);   MulMod(hh1, hh1, aa, HH);   conv(g1, gg1);   mul(g1, g1, p);   add(G, g, g1);   conv(h1, hh1);   mul(h1, h1, p);   add(H, h, h1);   Gout = G;   Hout = H;}staticvoid BuildTree(vec_long& link, vec_ZZX& v, vec_ZZX& w,               const vec_zz_pX& a){   long k = a.length();   if (k < 2) Error("bad arguments to BuildTree");   vec_zz_pX V, W;   V.SetLength(2*k-2);   W.SetLength(2*k-2);   link.SetLength(2*k-2);   long i, j, s;   long minp, mind;   for (i = 0; i < k; i++) {      V[i] = a[i];      link[i] = -(i+1);   }   for (j = 0; j < 2*k-4; j += 2) {      minp = j;      mind = deg(V[j]);      for (s = j+1; s < i; s++)         if (deg(V[s]) < mind) {            minp = s;            mind = deg(V[s]);         }      swap(V[j], V[minp]);      swap(link[j], link[minp]);      minp = j+1;      mind = deg(V[j+1]);      for (s = j+2; s < i; s++)         if (deg(V[s]) < mind) {            minp = s;            mind = deg(V[s]);         }      swap(V[j+1], V[minp]);      swap(link[j+1], link[minp]);      mul(V[i], V[j], V[j+1]);      link[i] = j;      i++;   }   zz_pX d;   for (j = 0; j < 2*k-2; j += 2) {      XGCD(d, W[j], W[j+1], V[j], V[j+1]);      if (!IsOne(d))         Error("relatively prime polynomials expected");   }   v.SetLength(2*k-2);   for (j = 0; j < 2*k-2; j++)      conv(v[j], V[j]);   w.SetLength(2*k-2);   for (j = 0; j < 2*k-2; j++)      conv(w[j], W[j]);}staticvoid RecTreeLift(const vec_long& link, vec_ZZX& v, vec_ZZX& w,                 const ZZ& p, const ZZX& f, long j, long inv){   if (j < 0) return;   if (inv)      HenselLift(v[j], v[j+1], w[j], w[j+1],                 f, v[j], v[j+1], w[j], w[j+1], p);   else      HenselLift1(v[j], v[j+1], f, v[j], v[j+1], w[j], w[j+1], p);   RecTreeLift(link, v, w, p, v[j], link[j], inv);   RecTreeLift(link, v, w, p, v[j+1], link[j+1], inv);}staticvoid TreeLift(const vec_long& link, vec_ZZX& v, vec_ZZX& w,               long e0, long e1, const ZZX& f, long inv)// lift from p^{e0} to p^{e1}{   ZZ p0, p1;   power(p0, zz_p::modulus(), e0);   power(p1, zz_p::modulus(), e1-e0);   ZZ_pBak bak;   bak.save();   ZZ_p::init(p1);   RecTreeLift(link, v, w, p0, f, v.length()-2, inv);   bak.restore();} void MultiLift(vec_ZZX& A, const vec_zz_pX& a, const ZZX& f, long e,               long verbose){   long k = a.length();   long i;   if (k < 2 || e < 1) Error("MultiLift: bad args");   if (!IsOne(LeadCoeff(f)))      Error("MultiLift: bad args");   for (i = 0; i < a.length(); i++)      if (!IsOne(LeadCoeff(a[i])))         Error("MultiLift: bad args");   if (e == 1) {      A.SetLength(k);      for (i = 0; i < k; i++)         conv(A[i], a[i]);      return;   }   vec_long E;   append(E, e);   while (e > 1) {      e = (e+1)/2;      append(E, e);   }   long l = E.length();   vec_ZZX v, w;   vec_long link;   double t;   if (verbose) {      cerr << "building tree...";      t = GetTime();   }   BuildTree(link, v, w, a);   if (verbose) cerr << (GetTime()-t) << "\n";   for (i = l-1; i > 0; i--) {      if (verbose) {         cerr << "lifting to " << E[i-1] << "...";         t = GetTime();      }            TreeLift(link, v, w, E[i], E[i-1], f, i != 1);      if (verbose) cerr << (GetTime()-t) << "\n";   }   A.SetLength(k);   for (i = 0; i < 2*k-2; i++) {      long t = link[i];      if (t < 0)         A[-(t+1)] = v[i];   }}staticvoid inplace_rev(ZZX& f){   long n = deg(f);   long i, j;   i = 0;   j = n;   while (i < j) {      swap(f.rep[i], f.rep[j]);      i++;      j--;   }   f.normalize();}long ZZXFac_InitNumPrimes = 7;long ZZXFac_MaxNumPrimes = 50;static void RecordPattern(vec_long& pat, vec_pair_zz_pX_long& fac){   long n = pat.length()-1;   long i;   for (i = 0; i <= n; i++)      pat[i] = 0;   long k = fac.length();   for (i = 0; i < k; i++) {      long d = fac[i].b;      long m = deg(fac[i].a)/d;      pat[d] = m;   }}staticlong NumFactors(const vec_long& pat){   long n = pat.length()-1;   long i;   long res = 0;   for (i = 0; i <= n; i++)       res += pat[i];   return res;}staticvoid CalcPossibleDegrees(ZZ& pd, const vec_long& pat){   long n = pat.length()-1;   set(pd);   long d, j;   ZZ t1;   for (d = 1; d <= n; d++)       for (j = 0; j < pat[d]; j++) {         LeftShift(t1, pd, d);         bit_or(pd, pd, t1);      }}static void CalcPossibleDegrees(vec_ZZ& S, const vec_ZZ_pX& fac, long k)// S[i] = possible degrees of the product of any subset of size k//        among fac[i...], encoded as a bit vector.      {   long r = fac.length();   S.SetLength(r);   if (r == 0)      return;   if (k < 1 || k > r)      Error("CalcPossibleDegrees: bad args");   long i, l;   ZZ old, t1;   set(S[r-1]);   LeftShift(S[r-1], S[r-1], deg(fac[r-1]));   for (i = r-2; i >= 0; i--) {      set(t1);      LeftShift(t1, t1, deg(fac[i]));      bit_or(S[i], t1, S[i+1]);   }   for (l = 2; l <= k; l++) {      old = S[r-l];      LeftShift(S[r-l], S[r-l+1], deg(fac[r-l]));      for (i = r-l-1; i >= 0; i--) {         LeftShift(t1, old, deg(fac[i]));         old = S[i];         bit_or(S[i], S[i+1], t1);      }   }}staticvec_zz_pX *SmallPrimeFactorization(LocalInfoT& LocalInfo, const ZZX& f,                            long verbose){   long n = deg(f);   long i;   double t;   LocalInfo.n = n;   long& NumPrimes = LocalInfo.NumPrimes;   NumPrimes = 0;   LocalInfo.NumFactors = 0;   // some sanity checking...   if (ZZXFac_InitNumPrimes < 1 || ZZXFac_InitNumPrimes > 10000)      Error("bad ZZXFac_InitNumPrimes");   if (ZZXFac_MaxNumPrimes < ZZXFac_InitNumPrimes || ZZXFac_MaxNumPrimes > 10000)      Error("bad ZZXFac_MaxNumPrimes");   LocalInfo.p.SetLength(ZZXFac_InitNumPrimes);   LocalInfo.pattern.SetLength(ZZXFac_InitNumPrimes);     // set bits 0..n of LocalInfo.PossibleDegrees    SetBit(LocalInfo.PossibleDegrees, n+1);   add(LocalInfo.PossibleDegrees, LocalInfo.PossibleDegrees, -1);   long minr = n+1;   long irred = 0;   vec_pair_zz_pX_long *bestfac = 0;   zz_pX *besth = 0;   vec_zz_pX *spfactors = 0;   zz_pContext bestp;   long bestp_index;   long maxroot = NextPowerOfTwo(deg(f))+1;   for (; NumPrimes < ZZXFac_InitNumPrimes;) {      long p = LocalInfo.s.next();      if (!p) Error("out of small primes");      if (divide(LeadCoeff(f), p)) {         if (verbose) cerr << "skipping " << p << "\n";         continue;      }      zz_p::init(p, maxroot);      zz_pX ff, ffp, d;      conv(ff, f);      MakeMonic(ff);      diff(ffp, ff);      GCD(d, ffp, ff);      if (!IsOne(d)) {         if (verbose)  cerr << "skipping " << p << "\n";         continue;      }      if (verbose) {         cerr << "factoring mod " << p << "...";         t = GetTime();      }      vec_pair_zz_pX_long thisfac;      zz_pX thish;       SFCanZass1(thisfac, thish, ff, 0);      LocalInfo.p[NumPrimes] = p;      vec_long& pattern = LocalInfo.pattern[NumPrimes];      pattern.SetLength(n+1);      RecordPattern(pattern, thisfac);      long r = NumFactors(pattern);            if (verbose) {         cerr << (GetTime()-t) << "\n";         cerr << "degree sequence: ";         for (i = 0; i <= n; i++)            if (pattern[i]) {               cerr << pattern[i] << "*" << i << " ";            }         cerr << "\n";      }      if (r == 1) {         irred = 1;         break;      }      // update admissibility info      ZZ pd;      CalcPossibleDegrees(pd, pattern);      bit_and(LocalInfo.PossibleDegrees, LocalInfo.PossibleDegrees, pd);      if (weight(LocalInfo.PossibleDegrees) == 2) {         irred = 1;         break;      }      if (r < minr) {         minr = r;         delete bestfac;         bestfac = NTL_NEW_OP vec_pair_zz_pX_long;         *bestfac = thisfac;         delete besth;         besth = NTL_NEW_OP zz_pX;         *besth = thish;         bestp.save();         bestp_index = NumPrimes;      }      NumPrimes++;   }