zzxfactoring.c

数值算法库for Unix

      t1 = acc;
      add(t1, t1, pb_eff_2);
      div(t1, t1, pb_eff);
      rem(t1, t1, pdelta);
      if (t1 > pdelta_2)
         sub(t1, t1, pdelta);

      C(i) = t1;
   }
}


static
void Compute_pb(vec_long& b,vec_ZZ& pb, long p, long d, 
                const ZZ& root_bound, long n)
{
   ZZ t1, t2;
   long i;

   t1 = 2*power(root_bound, d)*n;

   if (d == 1) {
      i = 0;
      t2 = 1;
   }
   else {
      i = b(d-1);
      t2 = pb(d-1);
   }

   while (t2 <= t1) {
      i++;
      t2 *= p;
   }

   b.SetLength(d);
   b(d) = i;

   pb.SetLength(d);
   pb(d) = t2;
}

static
void Compute_pdelta(long& delta, ZZ& pdelta, long p, long bit_delta)
{
   ZZ t1;
   long i;

   i = delta;
   t1 = pdelta;

   while (NumBits(t1) <= bit_delta) {
      i++;
      t1 *= p;
   }

   delta = i;
   pdelta = t1;
}

static
void BuildReductionMatrix(mat_ZZ& M, long& C, long r, long d, const ZZ& pdelta,
                          const vec_vec_ZZ& chop_vec, 
                          const mat_ZZ& B_L, long verbose)
{
   long s = B_L.NumRows();

   C = long( sqrt(double(d) * double(r)) / 2.0 ) + 1;

   M.SetDims(s+d, r+d);
   clear(M);


   long i, j, k;
   ZZ t1, t2;

   for (i = 1; i <= s; i++)
      for (j = 1; j <= r; j++)
         mul(M(i, j), B_L(i, j), C);

   ZZ pdelta_2;

   RightShift(pdelta_2, pdelta, 1);

   long maxbits = 0;

   for (i = 1; i <= s; i++)
      for (j = 1; j <= d; j++) {
         t1 = 0;
         for (k = 1; k <= r; k++) {
            mul(t2, B_L(i, k), chop_vec(k)(j));
            add(t1, t1, t2);
         }

         rem(t1, t1, pdelta);
         if (t1 > pdelta_2)
            sub(t1, t1, pdelta);

         maxbits = max(maxbits, NumBits(t1));

         M(i, j+r) = t1;
      }
  

   for (i = 1; i <= d; i++)
      M(i+s, i+r) = pdelta;

   if (verbose) 
      cerr << "ratio = " << double(maxbits)/double(NumBits(pdelta))
           << "; ";
}


static
void CutAway(mat_ZZ& B1, vec_ZZ& D, mat_ZZ& M, 
             long C, long r, long d)
{
   long k = M.NumRows();
   ZZ bnd = 4*to_ZZ(C)*to_ZZ(C)*to_ZZ(r) + to_ZZ(d)*to_ZZ(r)*to_ZZ(r);

   while (k >= 1 && 4*D[k] > bnd*D[k-1]) k--;

   mat_ZZ B2;

   B2.SetDims(k, r);
   long i, j;

   for (i = 1; i <= k; i++)
      for (j = 1; j <= r; j++)
         div(B2(i, j), M(i, j), C);

   M.kill(); // save space
   D.kill();

   ZZ det2;
   long rnk;

   rnk = image(det2, B2);

   B1.SetDims(rnk, r);
   for (i = 1; i <= rnk; i++)
      for (j = 1; j <= r; j++)
         B1(i, j) = B2(i + k - rnk, j);
}




static
long GotThem(vec_ZZX& factors, 
             const mat_ZZ& B_L,
             const vec_ZZ_pX& W, 
             const ZZX& f, 
             long bnd,
             long verbose)
{
   double tt0, tt1;
   ZZ det;
   mat_ZZ R;
   long s, r;
   long i, j, cnt;

   if (verbose) {
      cerr << "   checking A (s = " << B_L.NumRows() 
           << "): gauss...";
   }

   tt0 = GetTime();

   gauss(det, R, B_L);

   tt1 = GetTime();

   if (verbose) cerr << (tt1-tt0) << "; ";

   // check if condition A holds

   s = B_L.NumRows();
   r = B_L.NumCols();

   for (j = 0; j < r; j++) {
      cnt = 0;
      for (i = 0; i < s; i++) {
         if (R[i][j] == 0) continue;
         if (R[i][j] != det) {
            if (verbose) cerr << "failed.\n";
            return 0;
         }
         cnt++;
      }

      if (cnt != 1) {
         if (verbose) cerr << "failed.\n";
         return 0;
      }
   }

   if (verbose) {
      cerr << "passed.\n";
      cerr << "   checking B...";
   }

   // extract relevant information from R

   vec_vec_long I_vec;
   I_vec.SetLength(s);

   vec_long deg_vec;
   deg_vec.SetLength(s);

   for (i = 0; i < s; i++) {
      long dg = 0;

      for (j = 0; j < r; j++) {
         if (R[i][j] != 0) append(I_vec[i], j);
         dg += deg(W[j]);
      }

      deg_vec[i] = dg;
   }

   R.kill(); // save space


   // check if any candidate factor is the product of too few
   // modular factors

   for (i = 0; i < s; i++)
      if (I_vec[i].length() <= van_hoeij_card_thresh) {
         if (verbose) cerr << "X\n";
         return 0;
      }

   if (verbose) cerr << "1";


   // sort deg_vec, I_vec in order of increasing degree

   for (i = 0; i < s-1; i++)
      for (j = 0; j < s-1-i; j++)
         if (deg_vec[j] > deg_vec[j+1]) {
            swap(deg_vec[j], deg_vec[j+1]);
            swap(I_vec[j], I_vec[j+1]);
         }


   // perform constant term tests

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

   ZZ half_P;
   RightShift(half_P, ZZ_p::modulus(), 1);

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

   ZZ t1;

   for (i = 0; i < s; i++) {
      vec_long& I = I_vec[i];
      prod = lc;
      for (j = 0; j < I.length(); j++)
         mul(prod, prod, ConstTerm(W[I[j]]));

      t1 = rep(prod);
      if (t1 > half_P)
         sub(t1, t1, ZZ_p::modulus());

      if (!divide(ct, t1)) {
          if (verbose) cerr << "X\n";
          return 0;
      }
   }

   if (verbose) cerr << "2";


   // multiply out polynomials and perform size tests

   vec_ZZX fac;
   ZZ_pX gg;
   ZZX g;

   for (i = 0; i < s-1; i++) {
      vec_long& I = I_vec[i];
      mul(gg, W, I);
      mul(gg, gg, lc);
      BalCopy(g, gg);
      if (MaxBits(g) > bnd) {
         if (verbose) cerr << "X\n";
         return 0;
      }
      PrimitivePart(g, g);
      append(fac, g);
   }

   if (verbose) cerr << "3";


   // finally...trial division

   ZZX f1 = f;
   ZZX h;

   for (i = 0; i < s-1; i++) {
      if (!divide(h, f1, fac[i])) {
         cerr << "X\n";
         return 0;
      }

      f1 = h;
   }

   // got them!

   if (verbose) cerr << "$\n";

   append(factors, fac);
   append(factors, f1);

   return 1;
}


void AdditionalLifting(ZZ& P1, 
                       long& e1, 
                       vec_ZZX& w1, 
                       long p, 
                       long new_bound,
                       const ZZX& f, 
                       long doubling,
                       long verbose)
{
   long new_e1;

   if (doubling)
      new_e1 = max(2*e1, new_bound); // at least double e1
   else
      new_e1 = new_bound;

   if (verbose) {
      cerr << ">>> additional hensel lifting to " << new_e1 << "...\n";
   }

   ZZ new_P1;

   power(new_P1, p, new_e1);

   ZZX f1;
   ZZ t1, t2;
   long i;
   long n = deg(f);

   if (LeadCoeff(f) == 1)
      f1 = f;
   else if (LeadCoeff(f) == -1)
      negate(f1, f);
   else {
      rem(t1, LeadCoeff(f), new_P1);
      InvMod(t1, t1, new_P1);
      f1.rep.SetLength(n+1); 
      for (i = 0; i <= n; i++) {
         mul(t2, f.rep[i], t1);
         rem(f1.rep[i], t2, new_P1);
      }
   }

   zz_pBak bak;
   bak.save();

   zz_p::init(p, NextPowerOfTwo(n)+1);

   long r = w1.length();

   vec_zz_pX ww1;
   ww1.SetLength(r);
   for (i = 0; i < r; i++)
      conv(ww1[i], w1[i]);

   w1.kill();

   double tt0, tt1;

   tt0 = GetTime();

   MultiLift(w1, ww1, f1, new_e1, verbose);

   tt1 = GetTime();

   if (verbose) {
      cerr << "lifting time: " << (tt1-tt0) << "\n\n";
   }

   P1 = new_P1;
   e1 = new_e1;

   bak.restore();
}

static
void Compute_pb_eff(long& b_eff, ZZ& pb_eff, long p, long d, 
                    const ZZ& root_bound,  
                    long n, long ran_bits)
{
   ZZ t1, t2;
   long i;

   if (root_bound == 1)
      t1 = (to_ZZ(d)*to_ZZ(n)) << (ran_bits + 1);
   else
      t1 = (power(root_bound, d)*n) << (ran_bits + 2);

   i = 0;
   t2 = 1;

   while (t2 <= t1) {
      i++;
      t2 *= p;
   }

   b_eff = i;
   pb_eff = t2;
}



static
long d1_val(long bit_delta, long r, long s)
{
   return long( 0.30*double(r)*double(s)/double(bit_delta) ) + 1;
}




// Next comes van Hoeij's algorithm itself.
// Some notation that differs from van Hoeij's paper:
//   n = deg(f)
//   r = # modular factors
//   s = dim(B_L)  (gets smaller over time)
//   d = # traces used
//   d1 = number of "compressed" traces
//
// The algorithm starts with a "sparse" version of van Hoeij, so that
// at first the traces d = 1, 2, ... are used in conjunction with
// a d x d identity matrix for van Hoeij's matrix A.
// The number of "excess" bits used for each trace, bit_delta, is initially
// 2*r.
// 
// When d*bit_delta exceeds 0.25*r*s, we switch to 
// a "dense" mode, where we use only about 0.25*r*s "compressed" traces.
// These bounds follow from van Hoeij's heuristic estimates.
//
// In sparse mode, d and bit_delta increase exponentially (but gently).
// In dense mode, but d increases somewhat more aggressively,
// and bit_delta is increased more gently.


static
void FindTrueFactors_vH(vec_ZZX& factors, const ZZX& ff, 
                        const vec_ZZX& w, const ZZ& P, 
                        long p, long e,
                        LocalInfoT& LocalInfo,
                        long verbose,
                        long bnd)
{
   const long SkipSparse = 0;

   ZZ_pBak bak;
   bak.save();
   ZZ_p::init(P);

   long r = w.length();

   vec_ZZ_pX W;
   W.SetLength(r);

   long i, j;

   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() && 
      (k <= van_hoeij_card_thresh || W.length() <= van_hoeij_size_thresh)) {

      if (k <= 1)
         CardinalitySearch(factors, f, W, LocalInfo, k, bnd, verbose);
      else
         CardinalitySearch1(factors, f, W, LocalInfo, k, bnd, verbose);
      k++;
   }

   if (2*k > W.length()) {
      // rest is irreducible, so we're done

      append(factors, f);
   }
   else {

      // now we apply van Hoeij's algorithm proper to f
   
      double time_start, time_stop, lll_time, tt0, tt1;

      time_start = GetTime();
      lll_time = 0;
   
      if (verbose) {
         cerr << "\n\n*** starting knapsack procedure\n";
      }
   
      ZZ P1 = P;
      long e1 = e;    // invariant: P1 = p^{e1}
   
      r = W.length();
   
      vec_ZZX w1;
      w1.SetLength(r);
      for (i = 0; i < r; i++)
         conv(w1[i], W[i]);
   
      long n = deg(f);
   
      mat_ZZ B_L;            // van Hoeij's lattice
      ident(B_L, r);
   
      long d = 0;            // number of traces
      
      long bit_delta = 0;    // number of "excess" bits

      vec_long b;
      vec_ZZ pb;             // pb(i) = p^{b(i)}

      long delta = 0;
      ZZ pdelta = to_ZZ(1);  // pdelta = p^delta
      pdelta = 1;
   
      vec_vec_ZZ trace_vec;
      trace_vec.SetLength(r);
   
      vec_vec_ZZ chop_vec;
      chop_vec.SetLength(r);
   
      ZZ root_bound = RootBound(f);
   
      if (verbose) {
         cerr << "NumBits(root_bound) = " << NumBits(root_bound) << "\n";
      }

      long dense = 0;
      long ran_bits = 32;

      long loop_cnt = 0;


      long s = r;

      for (;;) {

         loop_cnt++;
   
         // if we are using the power hack, then we do not try too hard...
         // this is really a hack on a hack!

         if (ok_to_abandon && 
             ((d >= 2 && s > 128) || (d >= 3 && s > 32) || (d >= 4 && s > 8) ||
              d >= 5) ) {
            if (verbose) cerr << "   abandoning\n";
            append(factors, f);
            break;
         }

         long d_last, d_inc, d_index;

         d_last = d;

         // set d_inc: 

         if (!dense) {
            d_inc = 1 + d/8;
         }
         else {
            d_inc = 1 + d/4; 
         }

         d_inc = min(d_inc, n-1-d);
            
         d += d_inc;

         // set bit_delta:
   
         if (bit_delta == 0) {
            // set initial value...don't make it any smaller than 2*r

            bit_delta = 2*r; 
         }
         else {
            long extra_bits;

            if (!dense) {
               extra_bits = 1 + bit_delta/8;
            }
            else if (d_inc != 0) {
               if (d1_val(bit_delta, r, s) > 1)
                  extra_bits = 1 + bit_delta/16; 
               else
                  extra_bits = 0;
            }
            else
               extra_bits = 1 + bit_delta/8;

            bit_delta += extra_bits;
         }

         if (d > d1_val(bit_delta, r, s)) 
            dense = 1;
   
         Compute_pdelta(delta, pdelta, p, bit_delta);

         long d1;
         long b_eff;
         ZZ pb_eff;

         if (!dense) {
            for (d_in