#include <NTL/ZZ_pXFactoring.h>
#include <NTL/vec_ZZVec.h>
#include <NTL/fileio.h>
#include <NTL/FacVec.h>

#include <stdio.h>

#include <NTL/new.h>

NTL_START_IMPL





void SquareFreeDecomp(vec_pair_ZZ_pX_long& u, const ZZ_pX& ff)
{
   ZZ_pX f = ff;

   if (!IsOne(LeadCoeff(f)))
      Error("SquareFreeDecomp: bad args");

   ZZ_pX r, t, v, tmp1;
   long m, j, finished, done;

   u.SetLength(0);

   if (deg(f) == 0)
      return;

   m = 1;
   finished = 0;

   do {
      j = 1;
      diff(tmp1, f);
      GCD(r, f, tmp1);
      div(t, f, r);

      if (deg(t) > 0) {
         done = 0;
         do {
            GCD(v, r, t);
            div(tmp1, t, v);
            if (deg(tmp1) > 0) append(u, cons(tmp1, j*m));
            if (deg(v) > 0) {
               div(r, r, v);
               t = v;
               j++;
            }
            else
               done = 1;
         } while (!done);
         if (deg(r) == 0) finished = 1;
      }

      if (!finished) {
         /* r is a p-th power */
         long p, k, d;
         conv(p, ZZ_p::modulus());
         d = deg(r)/p;
         f.rep.SetLength(d+1);
         for (k = 0; k <= d; k++) 
            f.rep[k] = r.rep[k*p];
         m = m*p;
      }
   } while (!finished);
}
         


static
void NullSpace(long& r, vec_long& D, vec_ZZVec& M, long verbose)
{
   long k, l, n;
   long i, j;
   long pos;
   ZZ t1, t2;
   ZZ *x, *y;

   const ZZ& p = ZZ_p::modulus();

   n = M.length();

   D.SetLength(n);
   for (j = 0; j < n; j++) D[j] = -1;

   r = 0;

   l = 0;
   for (k = 0; k < n; k++) {

      if (verbose && k % 10 == 0) cerr << "+";

      pos = -1;
      for (i = l; i < n; i++) {
         rem(t1, M[i][k], p);
         M[i][k] = t1;
         if (pos == -1 && !IsZero(t1))
            pos = i;
      }

      if (pos != -1) {
         swap(M[pos], M[l]);

         // make M[l, k] == -1 mod p, and make row l reduced

         InvMod(t1, M[l][k], p);
         NegateMod(t1, t1, p);
         for (j = k+1; j < n; j++) {
            rem(t2, M[l][j], p);
            MulMod(M[l][j], t2, t1, p);
         }

         for (i = l+1; i < n; i++) {
            // M[i] = M[i] + M[l]*M[i,k]

            t1 = M[i][k];   // this is already reduced

            x = M[i].elts() + (k+1);
            y = M[l].elts() + (k+1);

            for (j = k+1; j < n; j++, x++, y++) {
               // *x = *x + (*y)*t1

               mul(t2, *y, t1);
               add(*x, *x, t2);
            }
         }

         D[k] = l;   // variable k is defined by row l
         l++;

      }
      else {
         r++;
      }
   }
}



static
void BuildMatrix(vec_ZZVec& M, long n, const ZZ_pX& g, const ZZ_pXModulus& F,
                 long verbose)
{
   long i, j, m;
   ZZ_pXMultiplier G;
   ZZ_pX h;

   ZZ t;
   sqr(t, ZZ_p::modulus());
   mul(t, t, n);

   long size = t.size();

   M.SetLength(n);
   for (i = 0; i < n; i++)
      M[i].SetSize(n, size);

   build(G, g, F);

   set(h);
   for (j = 0; j < n; j++) {
      if (verbose && j % 10 == 0) cerr << "+";

      m = deg(h);
      for (i = 0; i < n; i++) {
         if (i <= m)
            M[i][j] = rep(h.rep[i]);
         else
            clear(M[i][j]);
      }

      if (j < n-1)
         MulMod(h, h, G, F);
   }

   for (i = 0; i < n; i++)
      AddMod(M[i][i], M[i][i], -1, ZZ_p::modulus());

}



static
void RecFindRoots(vec_ZZ_p& x, const ZZ_pX& f)
{
   if (deg(f) == 0) return;

   if (deg(f) == 1) {
      long k = x.length();
      x.SetLength(k+1);
      negate(x[k], ConstTerm(f));
      return;
   }
      
   ZZ_pX h;

   ZZ_p r;
   ZZ p1;


   RightShift(p1, ZZ_p::modulus(), 1);
   
   {
      ZZ_pXModulus F;
      build(F, f);

      do {
         random(r);
         PowerXPlusAMod(h, r, p1, F);
         add(h, h, -1);
         GCD(h, h, f);
      } while (deg(h) <= 0 || deg(h) == deg(f));
   }

   RecFindRoots(x, h);
   div(h, f, h); 
   RecFindRoots(x, h);
}

void FindRoots(vec_ZZ_p& x, const ZZ_pX& ff)
{
   ZZ_pX f = ff;

   if (!IsOne(LeadCoeff(f)))
      Error("FindRoots: bad args");

   x.SetMaxLength(deg(f));
   x.SetLength(0);
   RecFindRoots(x, f);
}


static
void RandomBasisElt(ZZ_pX& g, const vec_long& D, const vec_ZZVec& M)
{
   ZZ t1, t2;

   long n = D.length();

   long i, j, s;

   g.rep.SetLength(n);

   vec_ZZ_p& v = g.rep;

   for (j = n-1; j >= 0; j--) {
      if (D[j] == -1)
         random(v[j]);
      else {
         i = D[j];

         // v[j] = sum_{s=j+1}^{n-1} v[s]*M[i,s]

         clear(t1);

         for (s = j+1; s < n; s++) {
            mul(t2, rep(v[s]), M[i][s]);
            add(t1, t1, t2);
         }

         conv(v[j], t1);
      }
   }
}



static
void split(ZZ_pX& f1, ZZ_pX& g1, ZZ_pX& f2, ZZ_pX& g2,
           const ZZ_pX& f, const ZZ_pX& g, 
           const vec_ZZ_p& roots, long lo, long mid)
{
   long r = mid-lo+1;

   ZZ_pXModulus F;
   build(F, f);

   vec_ZZ_p lroots(INIT_SIZE, r);
   long i;

   for (i = 0; i < r; i++)
      lroots[i] = roots[lo+i];


   ZZ_pX h, a, d;
   BuildFromRoots(h, lroots);
   CompMod(a, h, g, F);


   GCD(f1, a, f);
   
   div(f2, f, f1);

   rem(g1, g, f1);
   rem(g2, g, f2);
}

static
void RecFindFactors(vec_ZZ_pX& factors, const ZZ_pX& f, const ZZ_pX& g,
                    const vec_ZZ_p& roots, long lo, long hi)
{
   long r = hi-lo+1;

   if (r == 0) return;

   if (r == 1) {
      append(factors, f);
      return;
   }

   ZZ_pX f1, g1, f2, g2;

   long mid = (lo+hi)/2;

   split(f1, g1, f2, g2, f, g, roots, lo, mid);

   RecFindFactors(factors, f1, g1, roots, lo, mid);
   RecFindFactors(factors, f2, g2, roots, mid+1, hi);
}


static
void FindFactors(vec_ZZ_pX& factors, const ZZ_pX& f, const ZZ_pX& g,
                 const vec_ZZ_p& roots)
{
   long r = roots.length();

   factors.SetMaxLength(r);
   factors.SetLength(0);

   RecFindFactors(factors, f, g, roots, 0, r-1);
}

#if 0

static
void IterFindFactors(vec_ZZ_pX& factors, const ZZ_pX& f,
                     const ZZ_pX& g, const vec_ZZ_p& roots)
{
   long r = roots.length();
   long i;
   ZZ_pX h;

   factors.SetLength(r);

   for (i = 0; i < r; i++) {
      sub(h, g, roots[i]);
      GCD(factors[i], f, h);
   }
}

#endif


   

void SFBerlekamp(vec_ZZ_pX& factors, const ZZ_pX& ff, long verbose)
{
   ZZ_pX f = ff;

   if (!IsOne(LeadCoeff(f)))
      Error("SFBerlekamp: bad args");

   if (deg(f) == 0) {
      factors.SetLength(0);
      return;
   }

   if (deg(f) == 1) {
      factors.SetLength(1);
      factors[0] = f;
      return;
   }

   double t;

   const ZZ& p = ZZ_p::modulus();

   long n = deg(f);

   ZZ_pXModulus F;

   build(F, f);

   ZZ_pX g, h;

   if (verbose) { cerr << "computing X^p..."; t = GetTime(); }
   PowerXMod(g, p, F);
   if (verbose) { cerr << (GetTime()-t) << "\n"; }

   vec_long D;
   long r;

   vec_ZZVec M;

   if (verbose) { cerr << "building matrix..."; t = GetTime(); }
   BuildMatrix(M, n, g, F, verbose);
   if (verbose) { cerr << (GetTime()-t) << "\n"; }

   if (verbose) { cerr << "diagonalizing..."; t = GetTime(); }
   NullSpace(r, D, M, verbose);
   if (verbose) { cerr << (GetTime()-t) << "\n"; }


   if (verbose) cerr << "number of factors = " << r << "\n";

   if (r == 1) {
      factors.SetLength(1);
      factors[0] = f;
      return;
   }

   if (verbose) { cerr << "factor extraction..."; t = GetTime(); }

   vec_ZZ_p roots;

   RandomBasisElt(g, D, M);
   MinPolyMod(h, g, F, r);
   if (deg(h) == r) M.kill();
   FindRoots(roots, h);
   FindFactors(factors, f, g, roots);

   ZZ_pX g1;
   vec_ZZ_pX S, S1;
   long i;

   while (factors.length() < r) {
      if (verbose) cerr << "+";
      RandomBasisElt(g, D, M);
      S.kill();
      for (i = 0; i < factors.length(); i++) {
         const ZZ_pX& f = factors[i];
         if (deg(f) == 1) {
            append(S, f);
            continue;
         }
         build(F, f);
         rem(g1, g, F);
         if (deg(g1) <= 0) {
            append(S, f);
            continue;
         }
         MinPolyMod(h, g1, F, min(deg(f), r-factors.length()+1));
         FindRoots(roots, h);
         S1.kill();
         FindFactors(S1, f, g1, roots);
         append(S, S1);
      }
      swap(factors, S);
   }

   if (verbose) { cerr << (GetTime()-t) << "\n"; }

   if (verbose) {
      cerr << "degrees:";
      long i;
      for (i = 0; i < factors.length(); i++)
         cerr << " " << deg(factors[i]);
      cerr << "\n";
   }
}


void berlekamp(vec_pair_ZZ_pX_long& factors, const ZZ_pX& f, long verbose)
{
   double t;
   vec_pair_ZZ_pX_long sfd;
   vec_ZZ_pX x;

   if (!IsOne(LeadCoeff(f)))
      Error("berlekamp: bad args");

   
   if (verbose) { cerr << "square-free decomposition..."; t = GetTime(); }
   SquareFreeDecomp(sfd, f);
   if (verbose) cerr << (GetTime()-t) << "\n";

   factors.SetLength(0);

   long i, j;

   for (i = 0; i < sfd.length(); i++) {
      if (verbose) {
         cerr << "factoring multiplicity " << sfd[i].b 
              << ", deg = " << deg(sfd[i].a) << "\n";
      }

      SFBerlekamp(x, sfd[i].a, verbose);

      for (j = 0; j < x.length(); j++)
         append(factors, cons(x[j], sfd[i].b));
   }
}



static
void AddFactor(vec_pair_ZZ_pX_long& factors, const ZZ_pX& g, long d, long verbose)
{
   if (verbose)
      cerr << "degree=" << d << ", number=" << deg(g)/d << "\n";
   append(factors, cons(g, d));
}

static
void ProcessTable(ZZ_pX& f, vec_pair_ZZ_pX_long& factors, 
                  const ZZ_pXModulus& F, long limit, const vec_ZZ_pX& tbl,
                  long d, long verbose)

{
   if (limit == 0) return;

   if (verbose) cerr << "+";

   ZZ_pX t1;

   if (limit == 1) {
      GCD(t1, f, tbl[0]);
      if (deg(t1) > 0) {
         AddFactor(factors, t1, d, verbose);
         div(f, f, t1);
      }

      return;
   }

   long i;

   t1 = tbl[0];
   for (i = 1; i < limit; i++)
      MulMod(t1, t1, tbl[i], F);

   GCD(t1, f, t1);

   if (deg(t1) == 0) return;

   div(f, f, t1);

   ZZ_pX t2;

   i = 0;
   d = d - limit + 1;

   while (2*d <= deg(t1)) {
      GCD(t2, tbl[i], t1); 
      if (deg(t2) > 0) {
         AddFactor(factors, t2, d, verbose);
         div(t1, t1, t2);
      }

      i++;
      d++;
   }

   if (deg(t1) > 0)
      AddFactor(factors, t1, deg(t1), verbose);
}


void TraceMap(ZZ_pX& w, const ZZ_pX& a, long d, const ZZ_pXModulus& F, 
              const ZZ_pX& b)

{
   if (d < 0) Error("TraceMap: bad args");

   ZZ_pX y, z, t;

   z = b;
   y = a;
   clear(w);

   while (d) {
      if (d == 1) {
         if (IsZero(w)) 
            w = y;
         else {
            CompMod(w, w, z, F);
            add(w, w, y);
         }
      }
      else if ((d & 1) == 0) {
         Comp2Mod(z, t, z, y, z, F);
         add(y, t, y);
      }
      else if (IsZero(w)) {
         w = y;
         Comp2Mod(z, t, z, y, z, F);
         add(y, t, y);
      }
      else {
         Comp3Mod(z, t, w, z, y, w, z, F);
         add(w, w, y);
         add(y, t, y);
      }

      d = d >> 1;
   }
}


void PowerCompose(ZZ_pX& y, const ZZ_pX& h, long q, const ZZ_pXModulus& F)
{
   if (q < 0) Error("PowerCompose: bad args");

   ZZ_pX z(INIT_SIZE, F.n);
   long sw;

   z = h;
   SetX(y);

   while (q) {
      sw = 0;

      if (q > 1) sw = 2;
      if (q & 1) {
         if (IsX(y))
            y = z;
         else
            sw = sw | 1;
      }

      switch (sw) {
      case 0:
         break;

      case 1:
         CompMod(y, y, z, F);
         break;

      case 2:
         CompMod(z, z, z, F);
         break;

      case 3:
         Comp2Mod(y, z, y, z, z, F);
         break;
      }

      q = q >> 1;
   }
}


long ProbIrredTest(const ZZ_pX& f, long iter)