ratio[i] = to_ulong(aa);   }   InitTab(lookup_tab, ratio, r, k, thresh1, shamt_tab, pruning);   for (i = 0; i < r; i++)      for (j = 0; j < i; j++) {         mul(a, W[i].rep[deg(W[i])-1], W[j].rep[deg(W[j])-1]);         mul(a, a, lc);         LeftShift(aa, rep(a), NTL_BITS_PER_LONG);         div(aa, aa, p);         pair_ratio[i][j] = to_ulong(aa);      }   for (i = 0; i < r; i++) {      long m = deg(W[i]);      if (m >= 2) {         mul(a, W[i].rep[m-2], lc);         LeftShift(aa, rep(a), NTL_BITS_PER_LONG);         div(aa, aa, p);         pair_ratio[i][i] = to_ulong(aa);      }      else         pair_ratio[i][i] = 0;   }}static long SecondOrderTest(const vec_long& I_vec, const vec_vec_ulong& pair_ratio_vec,                     vec_ulong& sum_stack_vec, long& SumLen){   long k = I_vec.length();   const long *I = I_vec.elts();   unsigned long *sum_stack = sum_stack_vec.elts();   unsigned long sum, thresh1;   if (SumLen == 0) {      unsigned long epsilon = (1UL << (NTL_BITS_PER_LONG-ZZX_OVERLIFT));      unsigned long delta = (unsigned long) ((k*(k+1)) >> 1);      unsigned long thresh = epsilon + delta;      thresh1 = (epsilon << 1) + delta;      sum = thresh;      sum_stack[k] = thresh1;   }   else {      sum = sum_stack[SumLen-1];      thresh1 = sum_stack[k];   }   long i, j;   for (i = SumLen; i < k; i++) {      const unsigned long *p = pair_ratio_vec[I[i]].elts();      for (j = 0; j <= i; j++) {         sum += p[I[j]];      }      sum_stack[i] = sum;   }   SumLen = k-1;   return (sum <= thresh1);}staticZZ choose_fn(long r, long k){   ZZ a, b;   a = 1;    b = 1;   long i;   for (i = 0; i < k; i++) {      a *= r-i;      b *= k-i;   }   return a/b;}staticvoid PrintInfo(const char *s, const ZZ& a, const ZZ& b){   cerr << s << a << " / " << b << " = ";      double x = to_double(a)/to_double(b);   if (x == 0)       cerr << "0";    else {      int n;      double f;      f = frexp(x, &n);      cerr << f << "*2^" << n;   }   cerr << "\n";}staticvoid RemoveFactors1(vec_long& W, const vec_long& I, long r){   long k = I.length();   long i, j;   i = 0;   for (j = 0; j < r; j++) {      if (i < k && j == I[i])         i++;      else         swap(W[j-i], W[j]);    }}staticvoid RemoveFactors1(vec_vec_long& W, const vec_long& I, long r){   long k = I.length();   long i, j;   i = 0;   for (j = 0; j < r; j++) {      if (i < k && j == I[i])         i++;      else         swap(W[j-i], W[j]);    }   for (i = 0; i < r-k; i++)      RemoveFactors1(W[i], I, r);}// should this swap go in tools.h?// Maybe not...I don't want to pollute the interface too much more.inline void swap(unsigned long& a, unsigned long& b)     { unsigned long t;  t = a; a = b; b = t; }staticvoid RemoveFactors1(vec_ulong& W, const vec_long& I, long r){   long k = I.length();   long i, j;   i = 0;   for (j = 0; j < r; j++) {      if (i < k && j == I[i])         i++;      else         swap(W[j-i], W[j]);    }}staticvoid RemoveFactors1(vec_vec_ulong& W, const vec_long& I, long r){   long k = I.length();   long i, j;   i = 0;   for (j = 0; j < r; j++) {      if (i < k && j == I[i])         i++;      else         swap(W[j-i], W[j]);    }   for (i = 0; i < r-k; i++)      RemoveFactors1(W[i], I, r);}staticvoid RemoveFactors1(vec_ZZ_p& W, const vec_long& I, long r){   long k = I.length();   long i, j;   i = 0;   for (j = 0; j < r; j++) {      if (i < k && j == I[i])         i++;      else         swap(W[j-i], W[j]);   }}staticvoid SumCoeffs(ZZ& sum, const ZZX& a){   ZZ res;   res = 0;   long i;   long n = a.rep.length();   for (i = 0; i < n; i++)      res += a.rep[i];   sum = res;}staticvoid SumCoeffs(ZZ_p& sum, const ZZ_pX& a){   ZZ_p res;   res = 0;   long i;   long n = a.rep.length();   for (i = 0; i < n; i++)      res += a.rep[i];   sum = res;}staticlong ConstTermTest(const vec_ZZ_p& W,                   const vec_long& I,                  const ZZ& ct,                  const ZZ_p& lc,                  vec_ZZ_p& prod,                  long& ProdLen) {   long k = I.length();   ZZ_p t;   ZZ t1, t2;   long i;   if (ProdLen == 0) {      mul(prod[0], lc, W[I[0]]);      ProdLen++;   }   for (i = ProdLen; i < k; i++)      mul(prod[i], prod[i-1], W[I[i]]);   ProdLen = k-1;   // should make this a routine in ZZ_p   t1 = rep(prod[k-1]);   RightShift(t2, ZZ_p::modulus(), 1);   if (t1 > t2)      sub(t1, t1, ZZ_p::modulus());   return divide(ct, t1);}long ZZXFac_MaxPrune = 10;staticlong pruning_bnd(long r, long k){   double x = 0;    long i;   for (i = 0; i < k; i++) {      x += log(double(r-i)/double(k-i));   }   return long((x/log(2.0)) * 0.75);}staticlong shamt_tab_init(long pos, long card, long pruning, long thresh1_len){   double x = 1;   long i;   for (i = 0; i < card; i++) {      x *= double(pos-i)/double(card-i);   }   x *= pruning;  // this can be adjusted to control the density   if (pos <= 6) x *= 2;  // a little boost that costs very little         long t = long(ceil(log(x)/log(2.0)));   t = max(t, TBL_SHAMT);    t = min(t, NTL_BITS_PER_LONG-thresh1_len);   return NTL_BITS_PER_LONG-t;}// The following routine should only be called for k > 1,// and is only worth calling for k > 2.staticvoid CardinalitySearch1(vec_ZZX& factors, ZZX& f,                        vec_ZZ_pX& W,                        LocalInfoT& LocalInfo,                        long k,                       long bnd,                       long verbose){   double start_time, end_time;   if (verbose) {      start_time = GetTime();      cerr << "\n************ ";      cerr << "start cardinality " << k << "\n";   }   if (k <= 1) Error("internal error: call CardinalitySearch");   // This test is needed to ensure correcntes of "n-2" test   if (NumBits(k) > NTL_BITS_PER_LONG/2-2)      Error("Cardinality Search: k too large...");   vec_ZZ pdeg;   CalcPossibleDegrees(pdeg, W, k);   ZZ pd;   bit_and(pd, pdeg[0], LocalInfo.PossibleDegrees);   if (pd == 0) {      if (verbose) cerr << "skipping\n";      return;   }   vec_long I, D;   I.SetLength(k);   D.SetLength(k);   long r = W.length();   long initial_r = r;   vec_ulong ratio, ratio_sum;   ratio.SetLength(r);   ratio_sum.SetLength(k);   unsigned long epsilon = (1UL << (NTL_BITS_PER_LONG-ZZX_OVERLIFT));   unsigned long delta = (unsigned long) k;   unsigned long thresh = epsilon + delta;   unsigned long thresh1 = (epsilon << 1) + delta;   long thresh1_len = NumBits(long(thresh1));    long pruning;   pruning = min(r/2, ZZXFac_MaxPrune);   pruning = min(pruning, pruning_bnd(r, k));   pruning = min(pruning, NTL_BITS_PER_LONG-EXTRA_BITS-thresh1_len);   if (pruning <= 4) pruning = 0;   long init_pruning = pruning;   TBL_T ***lookup_tab = 0;   long **shamt_tab = 0;   if (pruning) {      typedef long *long_p;      long i, j;      shamt_tab = NTL_NEW_OP long_p[pruning+1];      if (!shamt_tab) Error("out of mem");      shamt_tab[0] = shamt_tab[1] = 0;      for (i = 2; i <= pruning; i++) {         long len = min(k-1, i);         shamt_tab[i] = NTL_NEW_OP long[len+1];         if (!shamt_tab[i]) Error("out of mem");         shamt_tab[i][0] = shamt_tab[i][1] = 0;         for (j = 2; j <= len; j++)            shamt_tab[i][j] = shamt_tab_init(i, j, pruning, thresh1_len);      }      typedef  TBL_T *TBL_T_p;      typedef  TBL_T **TBL_T_pp;      lookup_tab = NTL_NEW_OP TBL_T_pp[pruning+1];      if (!lookup_tab) Error("out of mem");      lookup_tab[0] = lookup_tab[1] = 0;      for (i = 2; i <= pruning; i++) {         long len = min(k-1, i);         lookup_tab[i] = NTL_NEW_OP TBL_T_p[len+1];         if (!lookup_tab[i]) Error("out of mem");         lookup_tab[i][0] = lookup_tab[i][1] = 0;         for (j = 2; j <= len; j++) {            lookup_tab[i][j] = NTL_NEW_OP TBL_T[((1L << (NTL_BITS_PER_LONG-shamt_tab[i][j]))+TBL_MSK) >> TBL_SHAMT];            if (!lookup_tab[i][j]) Error("out of mem");         }      }   }   if (verbose) {      cerr << "pruning = " << pruning << "\n";   }   vec_ZZ_p prod;   prod.SetLength(k);   long ProdLen;   vec_ZZ_p prod1;   prod1.SetLength(k);   long ProdLen1;   vec_ulong sum_stack;   sum_stack.SetLength(k+1);   long SumLen;   vec_long upd;   long i, state;   long cnt = 0;   ZZ ct;   mul(ct, ConstTerm(f), LeadCoeff(f));   ZZ_p lc;   conv(lc, LeadCoeff(f));   vec_vec_ulong pair_ratio;   pair_ratio.SetLength(r);   for (i = 0; i < r; i++)      pair_ratio[i].SetLength(r);   RatioInit1(ratio, W, lc, pruning, lookup_tab, pair_ratio, k, thresh1, shamt_tab);   ZZ c1;   SumCoeffs(c1, f);   mul(c1, c1, LeadCoeff(f));   vec_ZZ_p sum_coeffs;   sum_coeffs.SetLength(r);   for (i = 0; i < r; i++)      SumCoeffs(sum_coeffs[i], W[i]);   vec_long degv;   degv.SetLength(r);   for (i = 0; i < r; i++)      degv[i] = deg(W[i]);   ZZ_pX gg;   ZZX g, h;   I[0] = 0;     long loop_cnt = 0, degree_cnt = 0, n2_cnt = 0, sl_cnt = 0, ct_cnt = 0,         pl_cnt = 0, c1_cnt = 0, pl1_cnt = 0, td_cnt = 0;   ZZ loop_total, degree_total, n2_total, sl_total, ct_total,       pl_total, c1_total, pl1_total, td_total;   while (I[0] <= r-k) {      bit_and(pd, pdeg[I[0]], LocalInfo.PossibleDegrees);      if (IsZero(pd)) {         if (verbose) cerr << "skipping\n";         goto done;      }      unpack(upd, pd, LocalInfo.n);      D[0] = degv[I[0]];      ratio_sum[0] = ratio[I[0]] + thresh;      i = 1;      state = 0;      ProdLen = 0;      ProdLen1 = 0;      SumLen = 0;      for (;;) {         cnt++;         if (cnt > 2000000) {             if (verbose) {               loop_total += loop_cnt;  loop_cnt = 0;               degree_total += degree_cnt;  degree_cnt = 0;               n2_total += n2_cnt;  n2_cnt = 0;               sl_total += sl_cnt;  sl_cnt = 0;               ct_total += ct_cnt;  ct_cnt = 0;               pl_total += pl_cnt;  pl_cnt = 0;               c1_total += c1_cnt;  c1_cnt = 0;               pl1_total += pl1_cnt;  pl1_cnt = 0;               td_total += td_cnt;  td_cnt = 0;            }            cnt = 0;            UpdateLocalInfo(LocalInfo, pdeg, W, factors, f, k, verbose);            bit_and(pd, pdeg[I[0]], LocalInfo.PossibleDegrees);            if (IsZero(pd)) {               if (verbose) cerr << "skipping\n";               goto done;            }            unpack(upd, pd, LocalInfo.n);         }         if (i == k-1) {            unsigned long ratio_sum_last = ratio_sum[k-2];            long I_last = I[k-2];            {               long D_last = D[k-2];                  unsigned long rs;               long I_this;               long D_this;                  for (I_this = I_last+1; I_this < r; I_this++) {                  loop_cnt++;                     rs = ratio_sum_last + ratio[I_this];                  if (rs > thresh1) {                     cnt++;                     continue;                  }                  degree_cnt++;                     D_this = D_last + degv[I_this];                     if (!upd[D_this]) {                     cnt++;                     continue;                  }                     n2_cnt++;                  sl_cnt += (k-SumLen);                  I[k-1] = I_this;                  if (!SecondOrderTest(I, pair_ratio, sum_stack, SumLen)) {                     cnt += 2;                     continue;                  }                  c1_cnt++;                  pl1_cnt += (k-ProdLen1);                  if (!ConstTermTest(sum_coeffs, I, c1, lc, prod1, ProdLen1)) {                     cnt += 100;                     continue;                  }                  ct_cnt++;                  pl_cnt += (k-ProdLen);                  D[k-1] = D_this;                  if (!ConstTermTest(W, I, ct, lc, prod, ProdLen)) {                     cnt += 100;                     continue;                  }                  td_cnt++;                     if (verbose) {                     cerr << "+";                  }                     cnt += 1000;                     if (2*D[k-1] <= deg(f)) {                     mul(gg, W, I);                     mul(gg, gg, lc);                     BalCopy(g, gg);                     if(MaxBits(g) > bnd) {                        continue;                     }                     if (verbose) {                        cerr << "*";                     }                     PrimitivePart(g, g);                     if (!divide(h, f, g)) {                        continue;                     }                                       // factor found!                     append(factors, g);                     if (verbose) {                       cerr << "degree " << deg(g) << " factor found\n";                     }                     f = h;                     mul(ct, ConstTerm(f), LeadCoeff(f));                     conv(lc, LeadCoeff(f));                  }                  else {                     InvMul(gg, W, I);                     mul(gg, gg, lc);                     BalCopy(g, gg);                     if(MaxBits(g) > bnd) {                        continue;                     }                     if (verbose) {                        cerr << "*";                     }                     PrimitivePart(g, g);                     if (!divide(h, f, g)) {                        continue;                     }                           // factor found!                     append(factors, h);                     if (verbose) {                       cerr << "degree " << deg(h) << " factor found\n";                     }                     f = g;                     mul(ct, ConstTerm(f), LeadCoeff(f));                     conv(lc, LeadCoeff(f));                  }                        RemoveFactors(W, I);                  RemoveFactors1(degv, I, r);                  RemoveFactors1(sum_coeffs, I, r);                  RemoveFactors1(ratio, I, r);                  RemoveFactors1(pair_ratio, I, r);                  r = W.length();                  cnt = 0;                  pruning = min(pruning, r/2);                  if (pruning <= 4) pruning = 0;                  InitTab(lookup_tab, ratio, r, k, thresh1, shamt_tab, pruning);                  if (2*k > r)                      goto done;                  else                      goto restart;               } /* end of inner for loop */             }            i--;            state = 1;           }         else {            if (state == 0) {               long I_i = I[i-1] + 1;               I[i] = I_i;               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--;                  state = 1;               }               else {                  D[i] = D[i-1] + degv[I_i];                  ratio_sum[i] = ratio_sum[i-1] + ratio[I_i];                  i++;               }            }            else { // state == 1