zzxfactoring.cpp

数值算法库for Windows

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

      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)) {