#include <NTL/GF2XFactoring.h>

#include <NTL/new.h>

NTL_START_IMPL


long IterIrredTest(const GF2X& f)
{
   long df = deg(f);

   if (df <= 0) return 0;
   if (df == 1) return 1;

   GF2XModulus F;

   build(F, f);
   
   GF2X h;
   SetX(h);
   SqrMod(h, h, F);

   long i, d, limit, limit_sqr;
   GF2X g, X, t, prod;


   SetX(X);

   i = 0;
   g = h;
   d = 1;
   limit = 2;
   limit_sqr = limit*limit;

   set(prod);

   while (2*d <= df) {
      add(t, g, X);
      MulMod(prod, prod, t, F);
      i++;
      if (i == limit_sqr) {
         GCD(t, f, prod);
         if (!IsOne(t)) return 0;

         set(prod);
         limit++;
         limit_sqr = limit*limit;
         i = 0;
      }

      d = d + 1;
      if (2*d <= deg(f)) {
         SqrMod(g, g, F);
      }
   }

   if (i > 0) {
      GCD(t, f, prod);
      if (!IsOne(t)) return 0;
   }

   return 1;
}


void SquareFreeDecomp(vec_pair_GF2X_long& u, const GF2X& ff)
{
   GF2X f = ff;

   if (IsZero(f)) Error("SquareFreeDecomp: bad args");

   GF2X 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;
         p = 2;
         d = deg(r)/p;
         clear(f);
         for (k = 0; k <= d; k++)
            if (coeff(r, k*p) == 1)
               SetCoeff(f, k);

         m = m*p;
      }
   } while (!finished);
}
         



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

static
void ProcessTable(GF2X& f, vec_pair_GF2X_long& factors, 
                  const GF2XModulus& F, long limit, const vec_GF2X& tbl,
                  long d, long verbose)

{
   if (limit == 0) return;

   if (verbose) cerr << "+";

   GF2X 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);

   GF2X 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);
}


static
void TraceMap(GF2X& w, const GF2X& a, long d, const GF2XModulus& F)

{
   GF2X y, z;

   y = a;
   z = a;

   long i;

   for (i = 1; i < d; i++) {
      SqrMod(z, z, F);
      add(y, y, z);
   }

   w = y;
}


long GF2X_BlockingFactor = 40;

void DDF(vec_pair_GF2X_long& factors, const GF2X& ff, long verbose)
{
   GF2X f = ff;

   if (IsZero(f)) Error("DDF: bad args");

   factors.SetLength(0);

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

   if (deg(f) == 1) {
      AddFactor(factors, f, 1, verbose);
      return;
   }


   long GCDTableSize = GF2X_BlockingFactor;

   GF2XModulus F;
   build(F, f);

   long i, d, limit, old_n;
   GF2X g, X;


   vec_GF2X tbl(INIT_SIZE, GCDTableSize);

   SetX(X);

   i = 0;
   SqrMod(g, X, F);
   d = 1;
   limit = GCDTableSize;


   while (2*d <= deg(f)) {

      old_n = deg(f);
      add(tbl[i], g, X);
      i++;
      if (i == limit) {
         ProcessTable(f, factors, F, i, tbl, d, verbose);
         i = 0;
      }

      d = d + 1;
      if (2*d <= deg(f)) {
         // we need to go further

         if (deg(f) < old_n) {
            // f has changed 

            build(F, f);
            rem(g, g, F);
         }

         SqrMod(g, g, F);
      }
   }

   ProcessTable(f, factors, F, i, tbl, d-1, verbose);

   if (!IsOne(f)) AddFactor(factors, f, deg(f), verbose);
}




static
void EDFSplit(GF2X& f1, GF2X& f2, const GF2X& f, long d)
{
   GF2X a, g;
   GF2XModulus F;
   
   build(F, f);
   long n = F.n;

   do {
      random(a, n);
      TraceMap(g, a, d, F);
   } while (deg(g) <= 0);

   GCD(f1, f, g);
   div(f2, f, f1);
}

static
void RecEDF(vec_GF2X& factors, const GF2X& f, long d)
{
   if (deg(f) == d) {
      append(factors, f);
      return;
   }

   GF2X f1, f2;

   EDFSplit(f1, f2, f, d);
   RecEDF(factors, f1, d);
   RecEDF(factors, f2, d);
}
         

void EDF(vec_GF2X& factors, const GF2X& ff, long d, long verbose)

{
   GF2X f = ff;

   if (IsZero(f)) Error("EDF: bad args");

   long n = deg(f);
   long r = n/d;