zzxfactoring.c

密码大家Shoup写的数论算法c语言实现

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

   append(factors, f);

   bak.restore();
}




static
void ll_SFFactor(vec_ZZX& factors, const ZZX& ff, 
                 long verbose,
                 long bnd)

// input is primitive and square-free, with positive leading
// coefficient
{
   if (deg(ff) <= 1) {
      factors.SetLength(1);
      factors[0] = ff;
      if (verbose) {
         cerr << "*** SFFactor, trivial case 1.\n";
      }
      return;
   }

   // remove a factor of X, if necessary

   ZZX f;
   long xfac;
   long rev;

   double t;

   if (IsZero(ConstTerm(ff))) {
      RightShift(f, ff, 1);
      xfac = 1;
   }
   else {
      f = ff;
      xfac = 0;
   }

   // return a factor of X-1 if necessary

   long x1fac = 0;

   ZZ c1;
   SumCoeffs(c1, f);

   if (c1 == 0) {
      x1fac = 1;
      div(f, f, ZZX(1,1) - 1);
   }

   SumCoeffs(c1, f);

   if (deg(f) <= 1) {
      long r = 0;
      factors.SetLength(0);
      if (deg(f) > 0) {
         factors.SetLength(r+1);
         factors[r] = f;
         r++;
      }
      if (xfac) {
         factors.SetLength(r+1);
         SetX(factors[r]);
         r++;
      }

      if (x1fac) {
         factors.SetLength(r+1);
         factors[r] = ZZX(1,1) - 1;
         r++;
      }

      if (verbose) {
         cerr << "*** SFFactor: trivial case 2.\n";
      }

      return;
   }

   if (verbose) {
      cerr << "*** start SFFactor.\n";
   }

   // reverse f if this makes lead coefficient smaller

   ZZ t1, t2;

   abs(t1, LeadCoeff(f));
   abs(t2, ConstTerm(f));

   if (t1 > t2) {
      inplace_rev(f);
      rev = 1;
   }
   else 
      rev = 0;

   // obtain factorization modulo small primes

   if (verbose) {
      cerr << "factorization modulo small primes...\n";
      t = GetTime();
   }

   LocalInfoT LocalInfo;

   zz_pBak bak;
   bak.save();

   vec_zz_pX *spfactors =
       SmallPrimeFactorization(LocalInfo, f, verbose);

   if (!spfactors) {
      // f was found to be irreducible 

      bak.restore();

      if (verbose) {
         t = GetTime()-t;
         cerr << "small prime time: " << t << ", irreducible.\n";
      }

      if (rev)
         inplace_rev(f);

      long r = 0;

      factors.SetLength(r+1);
      factors[r] = f;
      r++;

      if (xfac) {
         factors.SetLength(r+1);
         SetX(factors[r]);
         r++;
      }

      if (x1fac) {
         factors.SetLength(r+1);
         factors[r] = ZZX(1,1) - 1;
         r++;
      }

      return;
   }

   if (verbose) {
      t = GetTime()-t;
      cerr << "small prime time: ";
      cerr << t << ", number of factors = " << spfactors->length() << "\n";
   }

   // prepare for Hensel lifting

   // first, calculate bit bound 

   long bnd1;
   long n = deg(f);
   long i;
   long e;
   ZZ P;
   long p;
   
   bnd1 = MaxBits(f) + (NumBits(n+1)+1)/2;

   if (!bnd || bnd1 < bnd)
      bnd = bnd1;

   i = n/2;
   while (!bit(LocalInfo.PossibleDegrees, i))
      i--;

   long lc_bnd = NumBits(LeadCoeff(f));

   long coeff_bnd = bnd + lc_bnd + i;

   long lift_bnd;

   lift_bnd = coeff_bnd + 15;  
   // +15 helps avoid trial divisions...can be any number >= 0

   lift_bnd = max(lift_bnd, bnd + lc_bnd + 2*NumBits(n) + ZZX_OVERLIFT);
   // facilitates "n-1" and "n-2" tests

   lift_bnd = max(lift_bnd, lc_bnd + NumBits(c1));
   // facilitates f(1) test

   lift_bnd += 2;
   // +2 needed to get inequalities right


   p = zz_p::modulus();

   e = long(double(lift_bnd)/(log(double(p))/log(double(2))));
   power(P, p, e);

   while (NumBits(P) <= lift_bnd) { 
      mul(P, P, p);
      e++;
   }

   if (verbose) {
      cerr << "lifting bound = " << lift_bnd << "\n";
      cerr << "Hensel lifting to exponent " << e << "...";
      t = GetTime();
   }

   // third, compute f1 so that it is monic and equal to f mod P

   ZZX f1;

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


   // Do Hensel lift

   vec_ZZX w;

   MultiLift(w, *spfactors, f1, e, verbose);


   if (verbose) {
      t = GetTime()-t;
      cerr << t << "\n";
   }

   // We're done with zz_p...restore

   delete spfactors;
   bak.restore();

   // search for true factors

   if (verbose) {
      cerr << "searching for true factors...\n";
      t = GetTime();
   }

   FindTrueFactors(factors, f, w, P, LocalInfo, 
                   verbose, coeff_bnd);

   if (verbose) {
      t = GetTime()-t;
      cerr << "factor search time " << t << "\n";
   }

   long r = factors.length();

   if (rev) {
      for (i = 0; i < r; i++) {
         inplace_rev(factors[i]);
         if (sign(LeadCoeff(factors[i])) < 0)
            negate(factors[i], factors[i]);
      }
   }

   if (xfac) {
      factors.SetLength(r+1);
      SetX(factors[r]);
      r++;
   }

   if (x1fac) {
      factors.SetLength(r+1);
      factors[r] = ZZX(1,1)-1;
      r++;
   }

   // that's it!!

   if (verbose) {
      cerr << "*** end SFFactor.  degree sequence:\n";
      for (i = 0; i < r; i++)
         cerr << deg(factors[i]) << " ";
      cerr << "\n";
   }
}



static 
long DeflationFactor(const ZZX& f)
{
   long n = deg(f);
   long m = 0;
   long i;

   for (i = 1; i <= n && m != 1; i++) {
      if (f.rep[i] != 0)
         m = GCD(m, i);
   }

   return m;
}

static
void inflate(ZZX& g, const ZZX& f, long m)
// input may not alias output
{
   long n = deg(f);
   long i;

   g = 0;
   for (i = n; i >= 0; i--) 
      SetCoeff(g, i*m, f.rep[i]);
}

static
void deflate(ZZX& g, const ZZX& f, long m)
// input may not alias output
{
   long n = deg(f);
   long i;

   g = 0;
   for (i = n; i >= 0; i -= m) 
      SetCoeff(g, i/m, f.rep[i]);
}

static
void MakeFacList(vec_long& v, long m)
{
   if (m <= 0) Error("internal error: MakeFacList");

   v.SetLength(0);

   long p = 2;
   while (m > 1) {
      while (m % p == 0)  {
         append(v, p);
         m = m / p;
      }

      p++;
   }
}

long ZZXFac_PowerHack = 1;

void SFFactor(vec_ZZX& factors, const ZZX& ff, 
              long verbose,
              long bnd)

// input is primitive and square-free, with positive leading
// coefficient

{
   if (ff == 0) 
      Error("SFFactor: bad args");

   if (deg(ff) <= 0) {
      factors.SetLength(0);
      return;
   }

   if (!ZZXFac_PowerHack) {
      ll_SFFactor(factors, ff, verbose, bnd);
      return;
   }

   long m = DeflationFactor(ff);

   if (m == 1) {
      if (verbose) {
         cerr << "SFFactor -- no deflation\n";
      }

      ll_SFFactor(factors, ff, verbose, bnd);
      return;
   }


   vec_long v;
   MakeFacList(v, m);
   long l = v.length();

   if (verbose) {
      cerr << "SFFactor -- deflation: " << v << "\n";
   }

   vec_ZZX res;
   res.SetLength(1);
   deflate(res[0], ff, m);

   long done;
   long j, k;

   done = 0;
   k = l-1;

   while (!done) {
      vec_ZZX res1;
      res1.SetLength(0);
      for (j = 0; j < res.length(); j++) {
         vec_ZZX res2;
         double t;
         if (verbose) {
            cerr << "begin - step " << k << ", " << j << "; deg = " 
                 << deg(res[j]) << "\n";
            t = GetTime();
         }

         ll_SFFactor(res2, res[j], verbose, k < 0 ? bnd : 0);

         if (verbose) {
            t = GetTime()-t;
            cerr << "end   - step " << k << ", " << j << "; time = "
                 << t << "\n\n";
         }

         append(res1, res2);
      }

      if (k < 0) {
         done = 1;
         swap(res, res1);
      }
      else {
         vec_ZZX res2;
         res2.SetLength(res1.length());
         for (j = 0; j < res1.length(); j++)
            inflate(res2[j], res1[j], v[k]);
         k--;
         swap(res, res2);
      }
   }

   factors = res;
}





void factor(ZZ& c,
            vec_pair_ZZX_long& factors,
            const ZZX& f,
            long verbose,
            long bnd)

{
   ZZX ff = f;

   if (deg(ff) <= 0) {
      c = ConstTerm(ff);
      factors.SetLength(0);
      return;
   }

   content(c, ff);
   divide(ff, ff, c);

   long bnd1 = MaxBits(ff) + (NumBits(deg(ff)+1)+1)/2;
   if (!bnd || bnd > bnd1)
      bnd = bnd1;

   vec_pair_ZZX_long sfd;

   double t;

   if (verbose) { cerr << "square-free decomposition..."; t = GetTime(); }
   SquareFreeDecomp(sfd, ff);
   if (verbose) cerr << (GetTime()-t) << "\n";

   factors.SetLength(0);

   vec_ZZX x;

   long i, j;

   for (i = 0; i < sfd.length(); i++) {
      if (verbose) {
         cerr << "factoring multiplicity " << sfd[i].b
              << ", deg = " << deg(sfd[i].a) << "\n";
         t = GetTime();
      }

      SFFactor(x, sfd[i].a, verbose, bnd);

      if (verbose) {
         t = GetTime()-t;
         cerr << "total time for multiplicity " 
              << sfd[i].b << ": " << t << "\n";
      }

      for (j = 0; j < x.length(); j++)
         append(factors, cons(x[j], sfd[i].b));
   }
}

NTL_END_IMPL