// h is replaced by the bivariate polynomial h*(X+Y) (mod f(X), g(Y)).{   long n = deg(g);   long k = h.length()-1;   if (k < 0) return;   if (k < n-1) {      h.SetLength(k+2);      h[k+1] = h[k];      for (long i = k; i >= 1; i--) {         MulByXMod(h[i], h[i], f);         add(h[i], h[i], h[i-1]);      }      MulByXMod(h[0], h[0], f);   }   else {      GF2EX b, t;      b = h[n-1];      for (long i = n-1; i >= 1; i--) {         mul(t, b, g.rep[i]);         MulByXMod(h[i], h[i], f);         add(h[i], h[i], h[i-1]);         add(h[i], h[i], t);      }      mul(t, b, g.rep[0]);      MulByXMod(h[0], h[0], f);      add(h[0], h[0], t);   }   // normalize   k = h.length()-1;   while (k >= 0 && IsZero(h[k])) k--;   h.SetLength(k+1);}staticvoid IrredCombine(GF2EX& x, const GF2EX& f, const GF2EX& g){   if (deg(f) < deg(g)) {      IrredCombine(x, g, f);      return;   }   // deg(f) >= deg(g)...not necessary, but maybe a little more   //                    time & space efficient   long df = deg(f);   long dg = deg(g);   long m = df*dg;   vec_GF2EX h(INIT_SIZE, dg);   long i;   for (i = 0; i < dg; i++) h[i].SetMaxLength(df);   h.SetLength(1);   set(h[0]);   vec_GF2E a;   a.SetLength(2*m);   for (i = 0; i < 2*m; i++) {      a[i] = ConstTerm(h[0]);      if (i < 2*m-1)         MulByXPlusY(h, f, g);   }   MinPolySeq(x, a, m);}staticvoid BuildPrimePowerIrred(GF2EX& f, long q, long e){   long n = power(q, e);   do {      random(f, n);      SetCoeff(f, n);   } while (!IterIrredTest(f));}staticvoid RecBuildIrred(GF2EX& f, long u, const FacVec& fvec){   if (fvec[u].link == -1)      BuildPrimePowerIrred(f, fvec[u].q, fvec[u].a);   else {      GF2EX g, h;      RecBuildIrred(g, fvec[u].link, fvec);      RecBuildIrred(h, fvec[u].link+1, fvec);      IrredCombine(f, g, h);   }}void BuildIrred(GF2EX& f, long n){   if (n <= 0)      Error("BuildIrred: n must be positive");   if (n >= (1L << (NTL_BITS_PER_LONG-4)))      Error("overflow in BuildIrred");   if (n == 1) {      SetX(f);      return;   }   FacVec fvec;   FactorInt(fvec, n);   RecBuildIrred(f, fvec.length()-1, fvec);}#if 0void BuildIrred(GF2EX& f, long n){   if (n <= 0)      Error("BuildIrred: n must be positive");   if (n >= (1L << (NTL_BITS_PER_LONG-4)))      Error("overflow in BuildIrred");   if (n == 1) {      SetX(f);      return;   }   GF2EX g;   do {      random(g, n);      SetCoeff(g, n);   } while (!IterIrredTest(g));   f = g;}#endifvoid BuildRandomIrred(GF2EX& f, const GF2EX& g){   GF2EXModulus G;   GF2EX h, ff;   build(G, g);   do {      random(h, deg(g));      IrredPolyMod(ff, h, G);   } while (deg(ff) < deg(g));   f = ff;}/************* NEW DDF ****************/long GF2EX_GCDTableSize = 4;char GF2EX_stem[256] = "";double GF2EXFileThresh = 256;static vec_GF2EX BabyStepFile;static vec_GF2EX GiantStepFile;static long use_files;staticdouble CalcTableSize(long n, long k){   double sz = GF2E::WordLength()+4;   sz = sz * sizeof(_ntl_ulong);   sz = sz * n;   sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_GF2E);   sz = sz * k;   sz = sz/1024;   return sz;}staticvoid GenerateBabySteps(GF2EX& h1, const GF2EX& f, const GF2EX& h, long k,                       long verbose){   double t;   if (verbose) { cerr << "generating baby steps..."; t = GetTime(); }   GF2EXModulus F;   build(F, f);   GF2EXArgument H;#if 0   double n2 = sqrt(double(F.n));   double n4 = sqrt(n2);   double n34 = n2*n4;   long sz = long(ceil(n34/sqrt(sqrt(2.0))));#else   long sz = 2*SqrRoot(F.n);#endif   build(H, h, F, sz);   h1 = h;   long i;   long HexOutput = GF2X::HexOutput;   GF2X::HexOutput = 1;   if (!use_files) {      BabyStepFile.kill();      BabyStepFile.SetLength(k-1);   }   for (i = 1; i <= k-1; i++) {      if (use_files) {         ofstream s;         OpenWrite(s, FileName(GF2EX_stem, "baby", i));         s << h1 << "\n";         s.close();      }      else         BabyStepFile(i) = h1;      CompMod(h1, h1, H, F);      if (verbose) cerr << "+";   }   if (verbose)      cerr << (GetTime()-t) << "\n";   GF2X::HexOutput = HexOutput;}staticvoid GenerateGiantSteps(const GF2EX& f, const GF2EX& h, long l, long verbose){   double t;   if (verbose) { cerr << "generating giant steps..."; t = GetTime(); }   GF2EXModulus F;   build(F, f);   GF2EXArgument H;#if 0   double n2 = sqrt(double(F.n));   double n4 = sqrt(n2);   double n34 = n2*n4;   long sz = long(ceil(n34/sqrt(sqrt(2.0))));#else   long sz = 2*SqrRoot(F.n);#endif   build(H, h, F, sz);   GF2EX h1;   h1 = h;   long i;   long HexOutput = GF2X::HexOutput;    GF2X::HexOutput = 1;   if (!use_files) {      GiantStepFile.kill();      GiantStepFile.SetLength(l);   }   for (i = 1; i <= l-1; i++) {      if (use_files) {         ofstream s;         OpenWrite(s, FileName(GF2EX_stem, "giant", i));         s << h1 << "\n";         s.close();      }      else         GiantStepFile(i) = h1;      CompMod(h1, h1, H, F);      if (verbose) cerr << "+";   }   if (use_files) {      ofstream s;      OpenWrite(s, FileName(GF2EX_stem, "giant", i));      s << h1 << "\n";      s.close();   }   else      GiantStepFile(i) = h1;   if (verbose)      cerr << (GetTime()-t) << "\n";   GF2X::HexOutput = HexOutput;}staticvoid FileCleanup(long k, long l){   if (use_files) {      long i;         for (i = 1; i <= k-1; i++)         remove(FileName(GF2EX_stem, "baby", i));         for (i = 1; i <= l; i++)         remove(FileName(GF2EX_stem, "giant", i));   }   else {      BabyStepFile.kill();      GiantStepFile.kill();   }}staticvoid NewAddFactor(vec_pair_GF2EX_long& u, const GF2EX& g, long m, long verbose){   long len = u.length();   u.SetLength(len+1);   u[len].a = g;   u[len].b = m;   if (verbose) {      cerr << "split " << m << " " << deg(g) << "\n";   }}   staticvoid NewProcessTable(vec_pair_GF2EX_long& u, GF2EX& f, const GF2EXModulus& F,                     vec_GF2EX& buf, long size, long StartInterval,                     long IntervalLength, long verbose){   if (size == 0) return;   GF2EX& g = buf[size-1];   long i;   for (i = 0; i < size-1; i++)      MulMod(g, g, buf[i], F);   GCD(g, f, g);   if (deg(g) == 0) return;   div(f, f, g);   long d = (StartInterval-1)*IntervalLength + 1;   i = 0;   long interval = StartInterval;   while (i < size-1 && 2*d <= deg(g)) {      GCD(buf[i], buf[i], g);      if (deg(buf[i]) > 0) {         NewAddFactor(u, buf[i], interval, verbose);         div(g, g, buf[i]);      }      i++;      interval++;      d += IntervalLength;   }   if (deg(g) > 0) {      if (i == size-1)         NewAddFactor(u, g, interval, verbose);      else         NewAddFactor(u, g, (deg(g)+IntervalLength-1)/IntervalLength, verbose);   }}staticvoid FetchGiantStep(GF2EX& g, long gs, const GF2EXModulus& F){   if (use_files) {      ifstream s;         OpenRead(s, FileName(GF2EX_stem, "giant", gs));         s >> g;      s.close();   }   else      g = GiantStepFile(gs);   rem(g, g, F);}staticvoid FetchBabySteps(vec_GF2EX& v, long k){   v.SetLength(k);   SetX(v[0]);   long i;   for (i = 1; i <= k-1; i++) {      if (use_files) {         ifstream s;         OpenRead(s, FileName(GF2EX_stem, "baby", i));         s >> v[i];         s.close();      }      else         v[i] = BabyStepFile(i);   }}      staticvoid GiantRefine(vec_pair_GF2EX_long& u, const GF2EX& ff, long k, long l,                 long verbose){   double t;   if (verbose) {      cerr << "giant refine...";      t = GetTime();   }   u.SetLength(0);   vec_GF2EX BabyStep;   FetchBabySteps(BabyStep, k);   vec_GF2EX buf(INIT_SIZE, GF2EX_GCDTableSize);   GF2EX f;   f = ff;   GF2EXModulus F;   build(F, f);   GF2EX g;   GF2EX h;   long size = 0;   long first_gs;   long d = 1;   while (2*d <= deg(f)) {      long old_n = deg(f);      long gs = (d+k-1)/k;      long bs = gs*k - d;      if (bs == k-1) {         size++;         if (size == 1) first_gs = gs;         FetchGiantStep(g, gs, F);         add(buf[size-1], g, BabyStep[bs]);      }      else {         add(h, g, BabyStep[bs]);         MulMod(buf[size-1], buf[size-1], h, F);      }      if (verbose && bs == 0) cerr << "+";      if (size == GF2EX_GCDTableSize && bs == 0) {         NewProcessTable(u, f, F, buf, size, first_gs, k, verbose);         if (verbose) cerr << "*";         size = 0;      }      d++;      if (2*d <= deg(f) && deg(f) < old_n) {         build(F, f);         long i;         for (i = 1; i <= k-1; i++)             rem(BabyStep[i], BabyStep[i], F);      }   }   if (size > 0) {      NewProcessTable(u, f, F, buf, size, first_gs, k, verbose);      if (verbose) cerr << "*";   }   if (deg(f) > 0)       NewAddFactor(u, f, 0, verbose);   if (verbose) {      t = GetTime()-t;      cerr << "giant refine time: " << t << "\n";   }}staticvoid IntervalRefine(vec_pair_GF2EX_long& factors, const GF2EX& ff,                    long k, long gs, const vec_GF2EX& BabyStep, long verbose){   vec_GF2EX buf(INIT_SIZE, GF2EX_GCDTableSize);   GF2EX f;   f = ff;   GF2EXModulus F;   build(F, f);   GF2EX g;   FetchGiantStep(g, gs, F);   long size = 0;   long first_d;   long d = (gs-1)*k + 1;   long bs = k-1;   while (bs >= 0 && 2*d <= deg(f)) {      long old_n = deg(f);      if (size == 0) first_d = d;      rem(buf[size], BabyStep[bs], F);      add(buf[size], buf[size], g);      size++;      if (size == GF2EX_GCDTableSize) {         NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose);         size = 0;      }      d++;      bs--;      if (bs >= 0 && 2*d <= deg(f) && deg(f) < old_n) {         build(F, f);         rem(g, g, F);      }   }   NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose);   if (deg(f) > 0)       NewAddFactor(factors, f, deg(f), verbose);}   staticvoid BabyRefine(vec_pair_GF2EX_long& factors, const vec_pair_GF2EX_long& u,                long k, long l, long verbose){   double t;   if (verbose) {      cerr << "baby refine...";      t = GetTime();   }   factors.SetLength(0);   vec_GF2EX BabyStep;   long i;   for (i = 0; i < u.length(); i++) {      const GF2EX& g = u[i].a;      long gs = u[i].b;      if (gs == 0 || 2*((gs-1)*k+1) > deg(g))         NewAddFactor(factors, g, deg(g), verbose);      else {         if (BabyStep.length() == 0)            FetchBabySteps(BabyStep, k);         IntervalRefine(factors, g, k, gs, BabyStep, verbose);      }   }   if (verbose) {      t = GetTime()-t;      cerr << "baby refine time: " << t << "\n";   }}                  void NewDDF(vec_pair_GF2EX_long& factors,            const GF2EX& f,            const GF2EX& h,            long verbose){   if (!IsOne(LeadCoeff(f)))      Error("NewDDF: bad args");   if (deg(f) == 0) {      factors.SetLength(0);      return;   }   if (deg(f) == 1) {      factors.SetLength(0);      append(factors, cons(f, 1));      return;   }   if (!GF2EX_stem[0])      sprintf(GF2EX_stem, "ddf-%ld", RandomBnd(10000));         long B = deg(f)/2;   long k = SqrRoot(B);   long l = (B+k-1)/k;   GF2EX h1;   if (CalcTableSize(deg(f), k + l - 1) > GF2EXFileThresh)      use_files = 1;   else      use_files = 0;   GenerateBabySteps(h1, f, h, k, verbose);   GenerateGiantSteps(f, h1, l, verbose);   vec_pair_GF2EX_long u;   GiantRefine(u, f, k, l, verbose);   BabyRefine(factors, u, k, l, verbose);   FileCleanup(k, l);}long IterComputeDegree(const GF2EX& h, const GF2EXModulus& F){   long n = deg(F);   if (n == 1 || IsX(h)) return 1;   long B = n/2;   long k = SqrRoot(B);   long l = (B+k-1)/k;   GF2EXArgument H;#if 0   double n2 = sqrt(double(n));   double n4 = sqrt(n2);   double n34 = n2*n4;   long sz = long(ceil(n34/sqrt(sqrt(2.0))));#else   long sz = 2*SqrRoot(F.n);#endif   build(H, h, F, sz);   GF2EX h1;   h1 = h;   vec_GF2EX baby;   baby.SetLength(k);   SetX(baby[0]);   long i;   for (i = 1; i <= k-1; i++) {      baby[i] = h1;      CompMod(h1, h1, H, F);      if (IsX(h1)) return i+1;   }   build(H, h1, F, sz);   long j;   for (j = 2; j <= l; j++) {      CompMod(h1, h1, H, F);      for (i = k-1; i >= 0; i--) {         if (h1 == baby[i])            return j*k-i;      }   }   return n;}NTL_END_IMPL