loop_cnt++;                     if (i < ProdLen)                  ProdLen = i;                     if (i < ProdLen1)                  ProdLen1 = i;                     if (i < SumLen)                  SumLen = i;               long I_i = (++I[i]);               if (i == 0) break;                  if (I_i > r-k+i) {                  i--;               }               else {                  long pruned;                  if (pruning && r-I_i <= pruning) {                     long pos = r-I_i;                     unsigned long rs = ratio_sum[i-1];                     unsigned long index1 = (rs >> shamt_tab[pos][k-i]);                     if (lookup_tab[pos][k-i][index1 >> TBL_SHAMT] & (1UL << (index1&TBL_MSK)))                        pruned = 0;                     else                        pruned = 1;                  }                  else                     pruned = 0;                      if (pruned) {                     i--;                  }                  else {                     D[i] = D[i-1] + degv[I_i];                     ratio_sum[i] = ratio_sum[i-1] + ratio[I_i];                     i++;                     state = 0;                  }               }            }         }      }      restart: ;   }   done:   if (lookup_tab) {      long i, j;      for (i = 2; i <= init_pruning; i++) {         long len = min(k-1, i);         for (j = 2; j <= len; j++) {            delete [] lookup_tab[i][j];         }         delete [] lookup_tab[i];      }      delete [] lookup_tab;   }   if (shamt_tab) {      long i;      for (i = 2; i <= init_pruning; i++) {         delete [] shamt_tab[i];      }      delete [] shamt_tab;   }   if (verbose) {       end_time = GetTime();      cerr << "\n************ ";      cerr << "end cardinality " << k << "\n";      cerr << "time: " << (end_time-start_time) << "\n";      ZZ loops_max = choose_fn(initial_r+1, k);      ZZ tuples_max = choose_fn(initial_r, k);      loop_total += loop_cnt;      degree_total += degree_cnt;      n2_total += n2_cnt;      sl_total += sl_cnt;      ct_total += ct_cnt;      pl_total += pl_cnt;      c1_total += c1_cnt;      pl1_total += pl1_cnt;      td_total += td_cnt;      cerr << "\n";      PrintInfo("loops: ", loop_total, loops_max);      PrintInfo("degree tests: ", degree_total, tuples_max);      PrintInfo("n-2 tests: ", n2_total, tuples_max);      cerr << "ave sum len: ";      if (n2_total == 0)          cerr << "--";      else         cerr << (to_double(sl_total)/to_double(n2_total));      cerr << "\n";      PrintInfo("f(1) tests: ", c1_total, tuples_max);      cerr << "ave prod len: ";      if (c1_total == 0)          cerr << "--";      else         cerr << (to_double(pl1_total)/to_double(c1_total));      cerr << "\n";      PrintInfo("f(0) tests: ", ct_total, tuples_max);      cerr << "ave prod len: ";      if (ct_total == 0)          cerr << "--";      else         cerr << (to_double(pl_total)/to_double(ct_total));      cerr << "\n";      PrintInfo("trial divs: ", td_total, tuples_max);   }}staticvoid FindTrueFactors(vec_ZZX& factors, const ZZX& ff,                      const vec_ZZX& w, const ZZ& P,                      LocalInfoT& LocalInfo,                     long verbose,                     long bnd){   ZZ_pBak bak;   bak.save();   ZZ_p::init(P);   long r = w.length();   vec_ZZ_pX W;   W.SetLength(r);   long i;   for (i = 0; i < r; i++)      conv(W[i], w[i]);   ZZX f;   f = ff;   long k;   k = 1;   factors.SetLength(0);   while (2*k <= W.length()) {      if (k <= 2)         CardinalitySearch(factors, f, W, LocalInfo, k, bnd, verbose);      else         CardinalitySearch1(factors, f, W, LocalInfo, k, bnd, verbose);      k++;   }   append(factors, f);   bak.restore();}staticvoid ll_SFFactor(vec_ZZX& factors, const ZZX& ff,                  long verbose,                 long bnd)// input is primitive and square-free, with positive leading// coefficient{   if (deg(ff) <= 1) {      factors.SetLength(1);      factors[0] = ff;      if (verbose) {         cerr << "*** SFFactor, trivial case 1.\n";      }      return;   }   // remove a factor of X, if necessary   ZZX f;   long xfac;   long rev;   double t;   if (IsZero(ConstTerm(ff))) {      RightShift(f, ff, 1);      xfac = 1;   }   else {      f = ff;      xfac = 0;   }   // return a factor of X-1 if necessary   long x1fac = 0;   ZZ c1;   SumCoeffs(c1, f);   if (c1 == 0) {      x1fac = 1;      div(f, f, ZZX(1,1) - 1);   }   SumCoeffs(c1, f);   if (deg(f) <= 1) {      long r = 0;      factors.SetLength(0);      if (deg(f) > 0) {         factors.SetLength(r+1);         factors[r] = f;         r++;      }      if (xfac) {         factors.SetLength(r+1);         SetX(factors[r]);         r++;      }      if (x1fac) {         factors.SetLength(r+1);         factors[r] = ZZX(1,1) - 1;         r++;      }      if (verbose) {         cerr << "*** SFFactor: trivial case 2.\n";      }      return;   }   if (verbose) {      cerr << "*** start SFFactor.\n";   }   // reverse f if this makes lead coefficient smaller   ZZ t1, t2;   abs(t1, LeadCoeff(f));   abs(t2, ConstTerm(f));   if (t1 > t2) {      inplace_rev(f);      rev = 1;   }   else       rev = 0;   // obtain factorization modulo small primes   if (verbose) {      cerr << "factorization modulo small primes...\n";      t = GetTime();   }   LocalInfoT LocalInfo;   zz_pBak bak;   bak.save();   vec_zz_pX *spfactors =       SmallPrimeFactorization(LocalInfo, f, verbose);   if (!spfactors) {      // f was found to be irreducible       bak.restore();      if (verbose) {         t = GetTime()-t;         cerr << "small prime time: " << t << ", irreducible.\n";      }      if (rev)         inplace_rev(f);      long r = 0;      factors.SetLength(r+1);      factors[r] = f;      r++;      if (xfac) {         factors.SetLength(r+1);         SetX(factors[r]);         r++;      }      if (x1fac) {         factors.SetLength(r+1);         factors[r] = ZZX(1,1) - 1;         r++;      }      return;   }   if (verbose) {      t = GetTime()-t;      cerr << "small prime time: ";      cerr << t << ", number of factors = " << spfactors->length() << "\n";   }   // prepare for Hensel lifting   // first, calculate bit bound    long bnd1;   long n = deg(f);   long i;   long e;   ZZ P;   long p;      bnd1 = MaxBits(f) + (NumBits(n+1)+1)/2;   if (!bnd || bnd1 < bnd)      bnd = bnd1;   i = n/2;   while (!bit(LocalInfo.PossibleDegrees, i))      i--;   long lc_bnd = NumBits(LeadCoeff(f));   long coeff_bnd = bnd + lc_bnd + i;   long lift_bnd;   lift_bnd = coeff_bnd + 15;     // +15 helps avoid trial divisions...can be any number >= 0   lift_bnd = max(lift_bnd, bnd + lc_bnd + 2*NumBits(n) + ZZX_OVERLIFT);   // facilitates "n-1" and "n-2" tests   lift_bnd = max(lift_bnd, lc_bnd + NumBits(c1));   // facilitates f(1) test   lift_bnd += 2;   // +2 needed to get inequalities right   p = zz_p::modulus();   e = long(double(lift_bnd)/(log(double(p))/log(double(2))));   power(P, p, e);   while (NumBits(P) <= lift_bnd) {       mul(P, P, p);      e++;   }   if (verbose) {      cerr << "lifting bound = " << lift_bnd << "\n";      cerr << "Hensel lifting to exponent " << e << "...";      t = GetTime();   }   // third, compute f1 so that it is monic and equal to f mod P   ZZX f1;   if (LeadCoeff(f) == 1)      f1 = f;   else if (LeadCoeff(f) == -1)      negate(f1, f);   else {      rem(t1, LeadCoeff(f), P);      if (sign(P) < 0)         Error("whoops!!!");      InvMod(t1, t1, P);      f1.rep.SetLength(n+1);      for (i = 0; i <= n; i++) {         mul(t2, f.rep[i], t1);         rem(f1.rep[i], t2, P);      }   }   // Do Hensel lift   vec_ZZX w;   MultiLift(w, *spfactors, f1, e, verbose);   if (verbose) {      t = GetTime()-t;      cerr << t << "\n";   }   // We're done with zz_p...restore   delete spfactors;   bak.restore();   // search for true factors   if (verbose) {      cerr << "searching for true factors...\n";      t = GetTime();   }   FindTrueFactors(factors, f, w, P, LocalInfo,                    verbose, coeff_bnd);   if (verbose) {      t = GetTime()-t;      cerr << "factor search time " << t << "\n";   }   long r = factors.length();   if (rev) {      for (i = 0; i < r; i++) {         inplace_rev(factors[i]);         if (sign(LeadCoeff(factors[i])) < 0)            negate(factors[i], factors[i]);      }   }   if (xfac) {      factors.SetLength(r+1);      SetX(factors[r]);      r++;   }   if (x1fac) {      factors.SetLength(r+1);      factors[r] = ZZX(1,1)-1;      r++;   }   // that's it!!   if (verbose) {      cerr << "*** end SFFactor.  degree sequence:\n";      for (i = 0; i < r; i++)         cerr << deg(factors[i]) << " ";      cerr << "\n";   }}static long DeflationFactor(const ZZX& f){   long n = deg(f);   long m = 0;   long i;   for (i = 1; i <= n && m != 1; i++) {      if (f.rep[i] != 0)         m = GCD(m, i);   }   return m;}staticvoid inflate(ZZX& g, const ZZX& f, long m)// input may not alias output{   long n = deg(f);   long i;   g = 0;   for (i = n; i >= 0; i--)       SetCoeff(g, i*m, f.rep[i]);}staticvoid deflate(ZZX& g, const ZZX& f, long m)// input may not alias output{   long n = deg(f);   long i;   g = 0;   for (i = n; i >= 0; i -= m)       SetCoeff(g, i/m, f.rep[i]);}staticvoid MakeFacList(vec_long& v, long m){   if (m <= 0) Error("internal error: MakeFacList");   v.SetLength(0);   long p = 2;   while (m > 1) {      while (m % p == 0)  {         append(v, p);         m = m / p;      }      p++;   }}long ZZXFac_PowerHack = 1;void SFFactor(vec_ZZX& factors, const ZZX& ff,               long verbose,              long bnd)// input is primitive and square-free, with positive leading// coefficient{   if (ff == 0)       Error("SFFactor: bad args");   if (deg(ff) <= 0) {      factors.SetLength(0);      return;   }   if (!ZZXFac_PowerHack) {      ll_SFFactor(factors, ff, verbose, bnd);      return;   }   long m = DeflationFactor(ff);   if (m == 1) {      if (verbose) {         cerr << "SFFactor -- no deflation\n";      }      ll_SFFactor(factors, ff, verbose, bnd);      return;   }   vec_long v;   MakeFacList(v, m);   long l = v.length();   if (verbose) {      cerr << "SFFactor -- deflation: " << v << "\n";   }   vec_ZZX res;   res.SetLength(1);   deflate(res[0], ff, m);   long done;   long j, k;   done = 0;   k = l-1;   while (!done) {      vec_ZZX res1;      res1.SetLength(0);      for (j = 0; j < res.length(); j++) {         vec_ZZX res2;         double t;         if (verbose) {            cerr << "begin - step " << k << ", " << j << "; deg = "                  << deg(res[j]) << "\n";            t = GetTime();         }         ll_SFFactor(res2, res[j], verbose, k < 0 ? bnd : 0);         if (verbose) {            t = GetTime()-t;            cerr << "end   - step " << k << ", " << j << "; time = "                 << t << "\n\n";         }         append(res1, res2);      }      if (k < 0) {         done = 1;         swap(res, res1);      }      else {         vec_ZZX res2;         res2.SetLength(res1.length());         for (j = 0; j < res1.length(); j++)            inflate(res2[j], res1[j], v[k]);         k--;         swap(res, res2);      }   }   factors = res;}void factor(ZZ& c,            vec_pair_ZZX_long& factors,            const ZZX& f,            long verbose,            long bnd){   ZZX ff = f;   if (deg(ff) <= 0) {      c = ConstTerm(ff);      factors.SetLength(0);      return;   }   content(c, ff);   divide(ff, ff, c);   long bnd1 = MaxBits(ff) + (NumBits(deg(ff)+1)+1)/2;   if (!bnd || bnd > bnd1)      bnd = bnd1;   vec_pair_ZZX_long sfd;   double t;   if (verbose) { cerr << "square-free decomposition..."; t = GetTime(); }   SquareFreeDecomp(sfd, ff);   if (verbose) cerr << (GetTime()-t) << "\n";   factors.SetLength(0);   vec_ZZX x;   long i, j;   for (i = 0; i < sfd.length(); i++) {      if (verbose) {         cerr << "factoring multiplicity " << sfd[i].b              << ", deg = " << deg(sfd[i].a) << "\n";         t = GetTime();      }      SFFactor(x, sfd[i].a, verbose, bnd);      if (verbose) {         t = GetTime()-t;         cerr << "total time for multiplicity "               << sfd[i].b << ": " << t << "\n";      }      for (j = 0; j < x.length(); j++)         append(factors, cons(x[j], sfd[i].b));   }}NTL_END_IMPL