GF2EX 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(GF2EX& y, const GF2EX& h, long q, const GF2EXModulus& F){   if (q < 0) Error("powerCompose: bad args");   GF2EX 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 GF2EX& f, long iter){   long n = deg(f);   if (n <= 0) return 0;   if (n == 1) return 1;   GF2EXModulus F;   build(F, f);   GF2EX b, r, s;   FrobeniusMap(b, F);   long all_zero = 1;   long i;   for (i = 0; i < iter; i++) {      random(r, n);      TraceMap(s, r, n, F, b);      all_zero = all_zero && IsZero(s);      if (deg(s) > 0) return 0;   }   if (!all_zero || (n & 1)) return 1;   PowerCompose(s, b, n/2, F);   return !IsX(s);}long GF2EX_BlockingFactor = 10;void DDF(vec_pair_GF2EX_long& factors, const GF2EX& ff, const GF2EX& hh,          long verbose){   GF2EX f = ff;   GF2EX h = hh;   if (!IsOne(LeadCoeff(f)))      Error("DDF: bad args");   factors.SetLength(0);   if (deg(f) == 0)      return;   if (deg(f) == 1) {      AddFactor(factors, f, 1, verbose);      return;   }   long CompTableSize = 2*SqrRoot(deg(f));    long GCDTableSize = GF2EX_BlockingFactor;   GF2EXModulus F;   build(F, f);   GF2EXArgument H;   build(H, h, F, min(CompTableSize, deg(f)));   long i, d, limit, old_n;   GF2EX g, X;   vec_GF2EX tbl(INIT_SIZE, GCDTableSize);   SetX(X);   i = 0;   g = h;   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(h, h, f);            rem(g, g, f);            build(H, h, F, min(CompTableSize, deg(f)));         }         CompMod(g, g, H, F);      }   }   ProcessTable(f, factors, F, i, tbl, d-1, verbose);   if (!IsOne(f)) AddFactor(factors, f, deg(f), verbose);}void RootEDF(vec_GF2EX& factors, const GF2EX& f, long verbose){   vec_GF2E roots;   double t;   if (verbose) { cerr << "finding roots..."; t = GetTime(); }   FindRoots(roots, f);   if (verbose) { cerr << (GetTime()-t) << "\n"; }   long r = roots.length();   factors.SetLength(r);   for (long j = 0; j < r; j++) {      SetX(factors[j]);      add(factors[j], factors[j], roots[j]);   }}staticvoid EDFSplit(vec_GF2EX& v, const GF2EX& f, const GF2EX& b, long d){   GF2EX a, g, h;   GF2EXModulus F;   vec_GF2E roots;      build(F, f);   long n = F.n;   long r = n/d;   random(a, n);   TraceMap(g, a, d, F, b);   MinPolyMod(h, g, F, r);   FindRoots(roots, h);   FindFactors(v, f, g, roots);}staticvoid RecEDF(vec_GF2EX& factors, const GF2EX& f, const GF2EX& b, long d,            long verbose){   vec_GF2EX v;   long i;   GF2EX bb;   if (verbose) cerr << "+";   EDFSplit(v, f, b, d);   for (i = 0; i < v.length(); i++) {      if (deg(v[i]) == d) {         append(factors, v[i]);      }      else {         GF2EX bb;         rem(bb, b, v[i]);         RecEDF(factors, v[i], bb, d, verbose);      }   }}         void EDF(vec_GF2EX& factors, const GF2EX& ff, const GF2EX& bb,         long d, long verbose){   GF2EX f = ff;   GF2EX b = bb;   if (!IsOne(LeadCoeff(f)))      Error("EDF: bad args");   long n = deg(f);   long r = n/d;   if (r == 0) {      factors.SetLength(0);      return;   }   if (r == 1) {      factors.SetLength(1);      factors[0] = f;      return;   }   if (d == 1) {      RootEDF(factors, f, verbose);      return;   }      double t;   if (verbose) {       cerr << "computing EDF(" << d << "," << r << ")...";       t = GetTime();    }   factors.SetLength(0);   RecEDF(factors, f, b, d, verbose);   if (verbose) cerr << (GetTime()-t) << "\n";}void SFCanZass(vec_GF2EX& factors, const GF2EX& ff, long verbose){   GF2EX f = ff;   if (!IsOne(LeadCoeff(f)))      Error("SFCanZass: bad args");   if (deg(f) == 0) {      factors.SetLength(0);      return;   }   if (deg(f) == 1) {      factors.SetLength(1);      factors[0] = f;      return;   }   factors.SetLength(0);   double t;      GF2EXModulus F;   build(F, f);   GF2EX h;   if (verbose) { cerr << "computing X^p..."; t = GetTime(); }   FrobeniusMap(h, F);   if (verbose) { cerr << (GetTime()-t) << "\n"; }   vec_pair_GF2EX_long u;   if (verbose) { cerr << "computing DDF..."; t = GetTime(); }   NewDDF(u, f, h, verbose);   if (verbose) {       t = GetTime()-t;       cerr << "DDF time: " << t << "\n";   }   GF2EX hh;   vec_GF2EX v;   long i;   for (i = 0; i < u.length(); i++) {      const GF2EX& g = u[i].a;      long d = u[i].b;      long r = deg(g)/d;      if (r == 1) {         // g is already irreducible         append(factors, g);      }      else {         // must perform EDF         if (d == 1) {            // root finding            RootEDF(v, g, verbose);            append(factors, v);         }         else {            // general case            rem(hh, h, g);            EDF(v, g, hh, d, verbose);            append(factors, v);         }      }   }}   void CanZass(vec_pair_GF2EX_long& factors, const GF2EX& f, long verbose){   if (!IsOne(LeadCoeff(f)))      Error("CanZass: bad args");   double t;   vec_pair_GF2EX_long sfd;   vec_GF2EX x;      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";      }      SFCanZass(x, sfd[i].a, verbose);      for (j = 0; j < x.length(); j++)         append(factors, cons(x[j], sfd[i].b));   }}void mul(GF2EX& f, const vec_pair_GF2EX_long& v){   long i, j, n;   n = 0;   for (i = 0; i < v.length(); i++)      n += v[i].b*deg(v[i].a);   GF2EX g(INIT_SIZE, n+1);   set(g);   for (i = 0; i < v.length(); i++)      for (j = 0; j < v[i].b; j++) {         mul(g, g, v[i].a);      }   f = g;}staticlong BaseCase(const GF2EX& h, long q, long a, const GF2EXModulus& F){   long b, e;   GF2EX lh(INIT_SIZE, F.n);   lh = h;   b = 1;   e = 0;   while (e < a-1 && !IsX(lh)) {      e++;      b *= q;      PowerCompose(lh, lh, q, F);   }   if (!IsX(lh)) b *= q;   return b;}staticvoid TandemPowerCompose(GF2EX& y1, GF2EX& y2, const GF2EX& h,                         long q1, long q2, const GF2EXModulus& F){   GF2EX z(INIT_SIZE, F.n);   long sw;   z = h;   SetX(y1);   SetX(y2);   while (q1 || q2) {      sw = 0;      if (q1 > 1 || q2 > 1) sw = 4;      if (q1 & 1) {         if (IsX(y1))            y1 = z;         else            sw = sw | 2;      }      if (q2 & 1) {         if (IsX(y2))            y2 = z;         else            sw = sw | 1;      }      switch (sw) {      case 0:         break;      case 1:         CompMod(y2, y2, z, F);         break;      case 2:         CompMod(y1, y1, z, F);         break;      case 3:         Comp2Mod(y1, y2, y1, y2, z, F);         break;      case 4:         CompMod(z, z, z, F);         break;      case 5:         Comp2Mod(z, y2, z, y2, z, F);         break;      case 6:         Comp2Mod(z, y1, z, y1, z, F);         break;      case 7:         Comp3Mod(z, y1, y2, z, y1, y2, z, F);         break;      }      q1 = q1 >> 1;      q2 = q2 >> 1;   }}staticlong RecComputeDegree(long u, const GF2EX& h, const GF2EXModulus& F,                      FacVec& fvec){   if (IsX(h)) return 1;   if (fvec[u].link == -1) return BaseCase(h, fvec[u].q, fvec[u].a, F);   GF2EX h1, h2;   long q1, q2, r1, r2;   q1 = fvec[fvec[u].link].val;    q2 = fvec[fvec[u].link+1].val;   TandemPowerCompose(h1, h2, h, q1, q2, F);   r1 = RecComputeDegree(fvec[u].link, h2, F, fvec);   r2 = RecComputeDegree(fvec[u].link+1, h1, F, fvec);   return r1*r2;}   long RecComputeDegree(const GF2EX& h, const GF2EXModulus& F)   // f = F.f is assumed to be an "equal degree" polynomial   // h = X^p mod f   // the common degree of the irreducible factors of f is computed{   if (F.n == 1 || IsX(h))       return 1;   FacVec fvec;   FactorInt(fvec, F.n);   return RecComputeDegree(fvec.length()-1, h, F, fvec);}void FindRoot(GF2E& root, const GF2EX& ff)// finds a root of ff.// assumes that ff is monic and splits into distinct linear factors{   GF2EXModulus F;   GF2EX h, h1, f;   GF2E r;   f = ff;      if (!IsOne(LeadCoeff(f)))      Error("FindRoot: bad args");   if (deg(f) == 0)      Error("FindRoot: bad args");   while (deg(f) > 1) {      build(F, f);      random(r);      clear(h);      SetCoeff(h, 1, r);      TraceMap(h, h, F);      GCD(h, h, f);      if (deg(h) > 0 && deg(h) < deg(f)) {         if (deg(h) > deg(f)/2)            div(f, f, h);         else            f = h;      }   }    root = ConstTerm(f);}staticlong power(long a, long e){   long i, res;   res = 1;   for (i = 1; i <= e; i++)      res = res * a;   return res;}staticlong IrredBaseCase(const GF2EX& h, long q, long a, const GF2EXModulus& F){   long e;   GF2EX X, s, d;   e = power(q, a-1);   PowerCompose(s, h, e, F);   SetX(X);   add(s, s, X);   GCD(d, F.f, s);   return IsOne(d);}staticlong RecIrredTest(long u, const GF2EX& h, const GF2EXModulus& F,                 const FacVec& fvec){   long  q1, q2;   GF2EX h1, h2;   if (IsX(h)) return 0;   if (fvec[u].link == -1) {      return IrredBaseCase(h, fvec[u].q, fvec[u].a, F);   }   q1 = fvec[fvec[u].link].val;    q2 = fvec[fvec[u].link+1].val;   TandemPowerCompose(h1, h2, h, q1, q2, F);   return RecIrredTest(fvec[u].link, h2, F, fvec)           && RecIrredTest(fvec[u].link+1, h1, F, fvec);}long DetIrredTest(const GF2EX& f){   if (deg(f) <= 0) return 0;   if (deg(f) == 1) return 1;   GF2EXModulus F;   build(F, f);      GF2EX h;   FrobeniusMap(h, F);   GF2EX s;   PowerCompose(s, h, F.n, F);   if (!IsX(s)) return 0;   FacVec fvec;   FactorInt(fvec, F.n);   return RecIrredTest(fvec.length()-1, h, F, fvec);}long IterIrredTest(const GF2EX& f){   if (deg(f) <= 0) return 0;   if (deg(f) == 1) return 1;   GF2EXModulus F;   build(F, f);      GF2EX h;   FrobeniusMap(h, F);   long CompTableSize = 2*SqrRoot(deg(f));   GF2EXArgument H;   build(H, h, F, CompTableSize);   long i, d, limit, limit_sqr;   GF2EX g, X, t, prod;   SetX(X);   i = 0;   g = h;   d = 1;   limit = 2;   limit_sqr = limit*limit;   set(prod);   while (2*d <= deg(f)) {      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)) {         CompMod(g, g, H, F);      }   }   if (i > 0) {      GCD(t, f, prod);      if (!IsOne(t)) return 0;   }   return 1;}staticvoid MulByXPlusY(vec_GF2EX& h, const GF2EX& f, const GF2EX& g)// h represents the bivariate polynomial h[0] + h[1]*Y + ... + h[n-1]*Y^k,// where the h[i]'s are polynomials in X, each of degree < deg(f),// and k < deg(g).