zzxfactoring.c

密码大家Shoup写的数论算法c语言实现

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


static
ZZ 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;
}

static
void 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";
}

static
void 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]); 
   }
}

static
void 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; }

static
void 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]); 
   }
}

static
void 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);
}


static
void 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]);
   }
}

static
void 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;
}

static
void 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;
}


static
long 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;



static
long 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);
}

static
long 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.


static
void 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