zzxfactoring.c

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

   if (!irred) {
      // delete best prime from LocalInfo
      swap(LocalInfo.pattern[bestp_index], LocalInfo.pattern[NumPrimes-1]);
      LocalInfo.p[bestp_index] = LocalInfo.p[NumPrimes-1];
      NumPrimes--;

      bestp.restore();

      spfactors = NTL_NEW_OP vec_zz_pX;

      if (verbose) {
         cerr << "p = " << zz_p::modulus() << ", completing factorization...";
         t = GetTime();
      }
      SFCanZass2(*spfactors, *bestfac, *besth, 0);
      if (verbose) {
         cerr << (GetTime()-t) << "\n";
      }
   }

   delete bestfac;
   delete besth;

   return spfactors;
}


static
long ConstTermTest(const vec_ZZ_pX& 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, ConstTerm(W[I[0]]));
      ProdLen++;
   }

   for (i = ProdLen; i < k; i++)
      mul(prod[i], prod[i-1], ConstTerm(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);
}

static
void BalCopy(ZZX& g, const ZZ_pX& G)
{
   const ZZ& p = ZZ_p::modulus();
   ZZ p2, t;
   RightShift(p2, p, 1);

   long n = G.rep.length();
   long i;

   g.rep.SetLength(n);
   for (i = 0; i < n; i++) {
      t = rep(G.rep[i]);
      if (t > p2) sub(t, t, p);
      g.rep[i] = t;
   }
}




static
void mul(ZZ_pX& g, const vec_ZZ_pX& W, const vec_long& I)
{
   vec_ZZ_pX w;
   long k = I.length();
   w.SetLength(k);
   long i;

   for (i = 0; i < k; i++)
      w[i] = W[I[i]];

   mul(g, w);
}




static
void InvMul(ZZ_pX& g, const vec_ZZ_pX& W, const vec_long& I)
{
   vec_ZZ_pX w;
   long k = I.length();
   long r = W.length();
   w.SetLength(r-k);
   long i, j;

   i = 0;
   for (j = 0; j < r; j++) {
      if (i < k && j == I[i])
         i++;
      else
         w[j-i] = W[j];
   } 

   mul(g, w);
}




static
void RemoveFactors(vec_ZZ_pX& W, const vec_long& I)
{
   long k = I.length();
   long r = W.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]); 
   }

   W.SetLength(r-k);
}

static
void unpack(vec_long& x, const ZZ& a, long n)
{
   x.SetLength(n+1);
   long i;

   for (i = 0; i <= n; i++)
      x[i] = bit(a, i);
}

static
void SubPattern(vec_long& p1, const vec_long& p2)
{
   long l = p1.length();

   if (p2.length() != l)
      Error("SubPattern: bad args");

   long i;

   for (i = 0; i < l; i++) {
      p1[i] -= p2[i];
      if (p1[i] < 0)
         Error("SubPattern: internal error");
   }
}

static
void UpdateLocalInfo(LocalInfoT& LocalInfo, vec_ZZ& pdeg,
                     const vec_ZZ_pX& W, const vec_ZZX& factors,
                     const ZZX& f, long k, long verbose)
{
   static long cnt = 0;

   if (verbose) {
      cnt = (cnt + 1) % 100;
      if (!cnt) cerr << "#";
   }

   double t;
   long i, j;

   if (LocalInfo.NumFactors < factors.length()) {
      zz_pBak bak;
      bak.save();

      vec_long pattern;
      pattern.SetLength(LocalInfo.n+1);

      ZZ pd;

      if (verbose) {
         cerr << "updating local info...";
         t = GetTime();
      }

      for (i = 0; i < LocalInfo.NumPrimes; i++) {
         zz_p::init(LocalInfo.p[i], NextPowerOfTwo(LocalInfo.n)+1);

         for (j = LocalInfo.NumFactors; j < factors.length(); j++) {
            vec_pair_zz_pX_long thisfac;
            zz_pX thish; 

            zz_pX ff;
            conv(ff, factors[j]);
            MakeMonic(ff);

            SFCanZass1(thisfac, thish, ff, 0);
            RecordPattern(pattern, thisfac);
            SubPattern(LocalInfo.pattern[i], pattern);
         }

         CalcPossibleDegrees(pd, LocalInfo.pattern[i]);
         bit_and(LocalInfo.PossibleDegrees, LocalInfo.PossibleDegrees, pd);

      }

      bak.restore();
      LocalInfo.NumFactors = factors.length();

      CalcPossibleDegrees(pdeg, W, k);

      if (verbose) cerr << (GetTime()-t) << "\n";
   }

   if (LocalInfo.NumPrimes < ZZXFac_MaxNumPrimes) {
      if (verbose)
         cerr << "adding a prime\n";

      zz_pBak bak;
      bak.save();

      for (;;) {
         long p = LocalInfo.s.next();
         if (!p)
            Error("UpdateLocalInfo: out of primes");

         if (divide(LeadCoeff(f), p)) {
            if (verbose) cerr << "skipping " << p << "\n";
            continue;
         }

         zz_p::init(p, NextPowerOfTwo(deg(f))+1);

         zz_pX ff, ffp, d;
   
         conv(ff, f);
         MakeMonic(ff);
         diff(ffp, ff);
   
         GCD(d, ffp, ff);
         if (!IsOne(d)) {
            if (verbose)  cerr << "skipping " << p << "\n";
            continue;
         }

         vec_pair_zz_pX_long thisfac;
         zz_pX thish;

         if (verbose) {
            cerr << "factoring mod " << p << "...";
            t = GetTime();
         }

         SFCanZass1(thisfac, thish, ff, 0);

         LocalInfo.p.SetLength(LocalInfo.NumPrimes+1);
         LocalInfo.pattern.SetLength(LocalInfo.NumPrimes+1);

         LocalInfo.p[LocalInfo.NumPrimes] = p;
         vec_long& pattern = LocalInfo.pattern[LocalInfo.NumPrimes];

         pattern.SetLength(LocalInfo.n+1);
         RecordPattern(pattern, thisfac);

         if (verbose) {
            cerr << (GetTime()-t) << "\n";
            cerr << "degree sequence: ";
            for (i = 0; i <= LocalInfo.n; i++)
               if (pattern[i]) {
                  cerr << pattern[i] << "*" << i << " ";
               }
            cerr << "\n";
         }

         ZZ pd;
         CalcPossibleDegrees(pd, pattern);
         bit_and(LocalInfo.PossibleDegrees, LocalInfo.PossibleDegrees, pd);

         LocalInfo.NumPrimes++;

         break;
      }

      bak.restore();
   }
}



const int ZZX_OVERLIFT = NTL_BITS_PER_LONG;
  // number of bits by which we "overlift"....this enables, in particular,
  // the "n-1" test.  
  // Must lie in the range 4..NTL_BITS_PER_LONG.


#define EXTRA_BITS (1)
// Any small number, like 1, 2 or 3, should be OK.


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

   vec_long I, D;
   I.SetLength(k);
   D.SetLength(k);

   long r = W.length();

   vec_ZZ_p prod;
   prod.SetLength(k);
   long ProdLen;

   vec_ZZ pdeg;
   CalcPossibleDegrees(pdeg, W, k);

   ZZ pd;
   vec_long upd;

   long i, state;

   long cnt = 0;

   ZZ ct;
   mul(ct, ConstTerm(f), LeadCoeff(f));

   ZZ_p lc;
   conv(lc, LeadCoeff(f));

   ZZ_pX gg;
   ZZX g, h;

   I[0] = 0;  

   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] = deg(W[I[0]]);
      i = 1;
      state = 0;
      ProdLen = 0;

      for (;;) {
         if (i < ProdLen)
            ProdLen = i;

         if (i == k) {
            // process indices I[0], ..., I[k-1]

            if (cnt > 2000000) { 
               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);
            }

            state = 1;  // default continuation state


            if (!upd[D[k-1]]) {
               i--;
               cnt++;
               continue;
            }

            if (!ConstTermTest(W, I, ct, lc, prod, ProdLen)) {
               i--;
               cnt += 100;
               continue;
            }

            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) {
                  i--;
                  continue;
               }
               if (verbose) {
                  cerr << "*";
               }
               PrimitivePart(g, g);
               if (!divide(h, f, g)) {
                  i--;
                  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) {
                  i--;
                  continue;
               }
               if (verbose) {
                  cerr << "*";
               }
               PrimitivePart(g, g);
               if (!divide(h, f, g)) {
                  i--;
                  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);
            r = W.length();
            cnt = 0;

            if (2*k > r) 
               goto done;
            else 
               break;
         }
         else if (state == 0) {
            I[i] = I[i-1] + 1;
            D[i] = D[i-1] + deg(W[I[i]]);
            i++;
         }
         else { // state == 1
            I[i]++;
            if (i == 0) break;

            if (I[i] > r-k+i)
               i--;
            else {
               D[i] = D[i-1] + deg(W[I[i]]);
               i++;
               state = 0;
            }
         }
      }
   }


   done: 


   if (verbose) {
      end_time = GetTime();
      cerr << "\n************ ";
      cerr << "end cardinality " << k << "\n";
      cerr << "time: " << (end_time-start_time) << "\n";
   }
}



typedef unsigned long TBL_T;

#if (NTL_BITS_PER_LONG >= 64)

// for 64-bit machines

#define TBL_MSK (63)
#define TBL_SHAMT (6)

#else

// for 32-bit machines

#define TBL_MSK (31)
#define TBL_SHAMT (5)

#endif


#if 0

// recursive version

static
void RecInitTab(TBL_T ***lookup_tab, long i, const vec_ulong& ratio, 
             long r, long k, unsigned long thresh1, long **shamt_tab,
             unsigned long sum, long card, long j)
{
   if (j >= i || card >= k-1) {
      if (card > 1) {
         long shamt = shamt_tab[i][card];
         unsigned long index1 = ((-sum) >> shamt);
         lookup_tab[i][card][index1 >> TBL_SHAMT] |= (1UL << (index1 & TBL_MSK));
         unsigned long index2 = ((-sum+thresh1) >> shamt);
         if (index1 != index2)
            lookup_tab[i][card][index2 >> TBL_SHAMT] |= (1UL << (index2 & TBL_MSK));

      }

      return;
   }


   RecInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab, sum, card, j+1);
   RecInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab, 
              sum+ratio[r-1-j], card+1, j+1);
}


static
void DoInitTab(TBL_T ***lookup_tab, long i, const vec_ulong& ratio, 
               long r, long k, unsigned long thresh1, long **shamt_tab)
{
   RecInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab, 0, 0, 0);
}

#else

// iterative version


static
void DoInitTab(TBL_T ***lookup_tab, long i, const vec_ulong& ratio,
               long r, long k, unsigned long thresh1, long **shamt_tab)
{
   vec_long sum_vec, card_vec, location_vec;
   sum_vec.SetLength(i+1);
   card_vec.SetLength(i+1);
   location_vec.SetLength(i+1);

   long j = 0;
   sum_vec[0] = 0;
   card_vec[0] = 0;

   unsigned long sum;
   long  card, location;

   location = 0;

   while (j >= 0) {
      sum = sum_vec[j];
      card = card_vec[j];

      switch (location) {

      case 0:

         if (j >= i || card >= k-1) {
            if (card > 1) {
               long shamt = shamt_tab[i][card];
               unsigned long index1 = ((-sum) >> shamt);
               lookup_tab[i][card][index1 >> TBL_SHAMT] |= (1UL << (index1 & TBL_MSK));
               unsigned long index2 = ((-sum+thresh1) >> shamt);
               if (index1 != index2)
                  lookup_tab[i][card][index2 >> TBL_SHAMT] |= (1UL << (index2 & TBL_MSK));
      
            }
      
            location = location_vec[j];
            j--;
            continue;
         }


         sum_vec[j+1] = sum;
         card_vec[j+1] = card;
         location_vec[j+1] = 1;
         j++;
         location = 0;
         continue;

      case 1:

         sum_vec[j+1] = sum+ratio[r-1-j];
         card_vec[j+1] = card+1;
         location_vec[j+1] = 2;
         j++;
         location = 0;
         continue;  

      case 2:

         location = location_vec[j];
         j--;
         continue;
      }
   }
}
         
#endif
   
   

static
void InitTab(TBL_T ***lookup_tab, const vec_ulong& ratio, long r, long k,
             unsigned long thresh1, long **shamt_tab, long pruning)
{
   long i, j, t;

   if (pruning) {
      for (i = 2; i <= pruning; i++) {
         long len = min(k-1, i);
         for (j = 2; j <= len; j++) {
            long ub = (((1L << (NTL_BITS_PER_LONG-shamt_tab[i][j])) 
                      + TBL_MSK) >> TBL_SHAMT); 
            for (t = 0; t < ub; t++)
               lookup_tab[i][j][t] = 0;
         }
   
         DoInitTab(lookup_tab, i, ratio, r, k, thresh1, shamt_tab);
      }
   }
}


static
void RatioInit1(vec_ulong& ratio, const vec_ZZ_pX& W, const ZZ_p& lc,
                long pruning, TBL_T ***lookup_tab, 
                vec_vec_ulong& pair_ratio, long k, unsigned long thresh1, 
                long **shamt_tab)
{
   long r = W.length();
   long i, j;

   ZZ_p a;

   ZZ p;
   p = ZZ_p::modulus();

   ZZ aa;

   for (i = 0; i < r; i++) {
      long m = deg(W[i]);
      mul(a, W[i].rep[m-1], lc);
      LeftShift(aa, rep(a), NTL_BITS_PER_LONG);
      div(aa, aa, p);