zzxfactoring.c

数值算法库for Unix

   
                  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

   if (d > n) {
      ZZ t1, t2;
      long i;

      t1 = 0;

      for (i = 1; i <= n; i++) {
         mul(t2, Tr(i + d - n - 1), f.rep[i-1]); 
         add(t1, t1, t2);
      }

      rem(t1, t1, P);
      NegateMod(t1, t1, P);
      Tr(d) = t1;
   }
   else {
      ZZ t1, t2;
      long i;

      mul(t1, f.rep[n-d], d);

      for (i = 1; i < d; i++) {
         mul(t2, Tr(i), f.rep[n-d+i]);
         add(t1, t1, t2);
      }

      rem(t1, t1, P);
      NegateMod(t1, t1, P);
      Tr(d) = t1;
   }
}

// Tr(1..d) are traces as computed above.
// C and pb have length at least d.
// For i = 1..d, pb(i) = p^{a_i} for a_i > 0.
// pdelta = p^delta for delta > 0.
// P = p^a for some a >= max{ a_i : i=1..d }.

// This routine computes C(1..d), where 
// C(i) = C_{a_i}^{a_i + delta}( Tr(i)*lc^i ) for i = 1..d.


void ChopTraces(vec_ZZ& C, const vec_ZZ& Tr, long d,
                const vec_ZZ& pb, const ZZ& pdelta, const ZZ& P, const ZZ& lc)
{
   if (d <= 0) Error("ChopTraces: internal error (1)");
   if (C.length() < d) Error("ChopTraces: internal error (2)");
   if (Tr.length() < d) Error("ChopTraces: internal error (3)");
   if (pb.length() < d) Error("ChopTraces: internal error (4)");
   if (P <= 1) Error("ChopTraces: internal error (5)");

   ZZ lcpow, lcred;
   lcpow = 1;
   rem(lcred, lc, P);

   ZZ pdelta_2;
   RightShift(pdelta_2, pdelta, 1);

   ZZ t1, t2;

   long i;
   for (i = 1; i <= d; i++) {
      MulMod(lcpow, lcpow, lcred, P);
      MulMod(t1, lcpow, Tr(i), P);

      RightShift(t2, pb(i), 1);
      add(t1, t1, t2);
      div(t1, t1, pb(i));
      rem(t1, t1, pdelta);
      if (t1 > pdelta_2)
         sub(t1, t1, pdelta);

      C(i) = t1;
   }
}


// Similar to above, but computes a linear combination of traces.


static
void DenseChopTraces(vec_ZZ& C, const vec_ZZ& Tr, long d, long d1, 
                     const ZZ& pb_eff, const ZZ& pdelta, const ZZ& P, 
                     const ZZ& lc, const mat_ZZ& A)
{

   ZZ pdelta_2;
   RightShift(pdelta_2, pdelta, 1);

   ZZ pb_eff_2;
   RightShift(pb_eff_2, pb_eff, 1);

   ZZ acc, t1, t2;

   long i, j;

   ZZ lcpow, lcred;
   rem(lcred, lc, P);

   for (i = 1; i <= d1; i++) {
      lcpow = 1;
      acc = 0;

      for (j = 1; j <= d; j++) {
         MulMod(lcpow, lcpow, lcred, P);
         MulMod(t1, lcpow, Tr(j), P);
         rem(t2, A(i, j), P);
         MulMod(t1, t1, t2, P);
         AddMod(acc, acc, t1, P);
      }