zzxfactoring.cpp

数值算法库for Windows

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



static
void 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 <= 1)
         CardinalitySearch(factors, f, W, LocalInfo, k, bnd, verbose);
      else
         CardinalitySearch1(factors, f, W, LocalInfo, k, bnd, verbose);
      k++;
   }

   append(factors, f);

   bak.restore();
}





/**********************************************************************\

                        van Hoeij's algorithm 

\**********************************************************************/



const long van_hoeij_size_thresh = 12; 
// Use van Hoeij's algorithm if number of modular factors exceeds this bound.
// Must be >= 1.

const long van_hoeij_card_thresh = 3;
// Switch to knapsack method if cardinality of candidate factors
// exceeds this bound.
// Must be >= 1.




// This routine assumes that the input f is a non-zero polynomial
// of degree n, and returns the value f(a).

static 
ZZ PolyEval(const ZZX& f, const ZZ& a)
{
   if (f == 0) Error("PolyEval: internal error");

   long n = deg(f);

   ZZ acc, t1, t2;
   long i;

   acc = f.rep[n];

   for (i = n-1; i >= 0; i--) {
      mul(t1, acc, a);
      add(acc, t1, f.rep[i]);
   }

   return acc;
}


// This routine assumes that the input f is a polynomial with non-zero constant
// term, of degree n, and with leading coefficient c; it returns 
// an upper bound on the absolute value of the roots of the
// monic, integer polynomial g(X) =  c^{n-1} f(X/c).

static 
ZZ RootBound(const ZZX& f)
{
   if (ConstTerm(f) == 0) Error("RootBound: internal error");

   long n = deg(f);

   ZZX g;
   long i;

   g = f;

   if (g.rep[n] < 0) negate(g.rep[n], g.rep[n]);
   for (i = 0; i < n; i++) {
      if (g.rep[i] > 0) negate(g.rep[i], g.rep[i]);
   }

   ZZ lb, ub, mb;


   lb = 0;

   ub = 1;
   while (PolyEval(g, ub) < 0) {
      ub = 2*ub;
   }

   // lb < root <= ub

   while (ub - lb > 1) {
      ZZ mb = (ub + lb)/2;

      if (PolyEval(g, mb) < 0) 
         lb = mb;
      else 
         ub = mb;
   }

   return ub*g.rep[n];
}


// This routine takes as input an n x m integer matrix M, where the rows of M 
// are assumed to be linearly independent.
// It is also required that both n and m are non-zero.
// It computes an integer d, along with an n x m matrix R, such that
// R*d^{-1} is the reduced row echelon form of M.
// The routine is probabilistic: the output is always correct, but the
// routine may abort the program with negligible probability
// (specifically, if GenPrime returns a composite, and the modular
// gauss routine can't invert a non-zero element).

static
void gauss(ZZ& d_out, mat_ZZ& R_out, const mat_ZZ& M)
{
   long n = M.NumRows();
   long m = M.NumCols();

   if (n == 0 || m == 0) Error("gauss: internal error");

   zz_pBak bak;
   bak.save();

   for (;;) {
      long p = GenPrime_long(NTL_SP_NBITS);
      zz_p::init(p);

      mat_zz_p MM;
      conv(MM, M);

      long r = gauss(MM);
      if (r < n) continue;

      // compute pos(1..n), so that pos(i) is the index 
      // of the i-th pivot column

      vec_long pos;
      pos.SetLength(n);

      long i, j;
      for (i = j = 1; i <= n; i++) {
         while (MM(i, j) == 0) j++;
         pos(i) = j;
         j++;
      } 

      // compute the n x n sub-matrix consisting of the
      // pivot columns of M

      mat_ZZ S;
      S.SetDims(n, n);

      for (i = 1; i <= n; i++)
         for (j = 1; j <= n; j++)
            S(i, j) = M(i, pos(j));

      mat_ZZ S_inv;
      ZZ d;

      inv(d, S_inv, S);
      if (d == 0) continue;

      mat_ZZ R;
      mul(R, S_inv, M);

      // now check that R is of the right form, which it will be
      // if we were not unlucky

      long OK = 1;

      for (i = 1; i <= n && OK; i++) {
         for (j = 1; j < pos(i) && OK; j++)
            if (R(i, j) != 0) OK = 0;

         if (R(i, pos(i)) != d) OK = 0;

         for (j = 1; j < i && OK; j++)
            if (R(j, pos(i)) != 0) OK = 0;
      }

      if (!OK) continue;

      d_out = d;
      R_out = R;
      break;
   }
}


// The input polynomial f should be monic, and deg(f) > 0.
// The input P should be > 1.
// Tr.length() >= d, and Tr(i), for i = 1..d-1, should be the
// Tr_i(f) mod P (in van Hoeij's notation).
// The quantity Tr_d(f) mod P is computed, and stored in Tr(d).


void ComputeTrace(vec_ZZ& Tr, const ZZX& f, long d, const ZZ& P)
{
   long n = deg(f);

   // check arguments

   if (n <= 0 || LeadCoeff(f) != 1) 
      Error("ComputeTrace: internal error (1)");

   if (d <= 0)
      Error("ComputeTrace: internal error (2)");

   if (Tr.length() < d)
      Error("ComputeTrace: internal error (3)");

   if (P <= 1)
      Error("ComputeTrace: internal error (4)");

   // treat d > deg(f) separately