zzxfactoring.cpp

数值算法库for Windows

         if (verbose) cerr << "skipping " << p << "\n";
         continue;
      }
      zz_p::init(p, maxroot);

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


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

      vec_pair_zz_pX_long thisfac;
      zz_pX thish; 

      SFCanZass1(thisfac, thish, ff, 0);

      LocalInfo.p[NumPrimes] = p;

      vec_long& pattern = LocalInfo.pattern[NumPrimes];
      pattern.SetLength(n+1);

      RecordPattern(pattern, thisfac);
      long r = NumFactors(pattern);
      
      if (verbose) {
         cerr << (GetTime()-t) << "\n";
         cerr << "degree sequence: ";
         for (i = 0; i <= n; i++)
            if (pattern[i]) {
               cerr << pattern[i] << "*" << i << " ";
            }
         cerr << "\n";
      }

      if (r == 1) {
         irred = 1;
         break;
      }

      // update admissibility info

      ZZ pd;

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

      if (weight(LocalInfo.PossibleDegrees) == 2) {
         irred = 1;
         break;
      }


      if (r < minr) {
         minr = r;
         delete bestfac;
         bestfac = NTL_NEW_OP vec_pair_zz_pX_long;
         *bestfac = thisfac;
         delete besth;
         besth = NTL_NEW_OP zz_pX;
         *besth = thish;
         bestp.save();
         bestp_index = NumPrimes;
      }

      NumPrimes++;
   }

   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 (!ZZXFac_van_Hoeij && LocalInfo.NumPrimes + 1 < 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