zzxfactoring.c

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

#include <NTL/ZZXFactoring.h>
#include <NTL/lzz_pXFactoring.h>
#include <NTL/vec_vec_long.h>
#include <NTL/vec_vec_ulong.h>
#include <NTL/vec_double.h>

#include <NTL/new.h>

NTL_START_IMPL


struct LocalInfoT {
   long n;
   long NumPrimes;
   long NumFactors;
   vec_long p;
   vec_vec_long pattern;
   ZZ PossibleDegrees;
   PrimeSeq s;
};



static
void mul(ZZ_pX& x, vec_ZZ_pX& a)
// this performs multiplications in close-to-optimal order,
// and kills a in the process
{
   long n = a.length();

   // first, deal with some trivial cases

   if (n == 0) {
      set(x);
      a.kill();
      return;
   }
   else if (n == 1) {
      x = a[0];
      a.kill();
      return;
   }

   long i, j;

   // assume n > 1 and all a[i]'s are nonzero

   // sort into non-increasing degrees

   for (i = 1; i <= n - 1; i++)
      for (j = 0; j <= n - i - 1; j++)
         if (deg(a[j]) < deg(a[j+1]))
            swap(a[j], a[j+1]);

   ZZ_pX g;

   while (n > 1) {
      // replace smallest two poly's by their product
      mul(g, a[n-2], a[n-1]);
      a[n-2].kill();
      a[n-1].kill();
      swap(g, a[n-2]);
      n--;

      // re-establish order

      i = n-1;
      while (i > 0 && deg(a[i-1]) < deg(a[i])) {
         swap(a[i-1], a[i]);
         i--;
      }
   }

   x = a[0];

   a[0].kill();
   a.SetLength(0);
}


void mul(ZZX& x, const vec_pair_ZZX_long& a)
{
   long l = a.length();
   ZZX res;
   long i, j;

   set(res);
   for (i = 0; i < l; i++)
      for (j = 0; j < a[i].b; j++)
         mul(res, res, a[i].a);

   x = res;
}


void SquareFreeDecomp(vec_pair_ZZX_long& u, const ZZX& ff)
// input is primitive 
{
   ZZX f = ff;

   ZZX d, v, w, s, t1;
   long i;

   u.SetLength(0);

   if (deg(f) <= 0)
      return;

   diff(t1, f);
   GCD(d, f, t1);

   if (deg(d) == 0) {
      append(u, cons(f, 1));
      return;
   }

   divide(v, f, d); 
   divide(w, t1, d);
   i = 0;

   for (;;) {
      i = i + 1;

      diff(t1, v);
      sub(s, w, t1);

      if (IsZero(s)) {
         if (deg(v) != 0) append(u, cons(v, i));
         return;
      }

      GCD(d, v, s);
      divide(v, v, d);
      divide(w, s, d);

      if (deg(d) != 0) append(u, cons(d, i));
   }
}




static
void HenselLift(ZZX& Gout, ZZX& Hout, ZZX& Aout, ZZX& Bout,
                const ZZX& f, const ZZX& g, const ZZX& h,
                const ZZX& a, const ZZX& b, const ZZ& p) 
{
   ZZX c, g1, h1, G, H, A, B;

   mul(c, g, h);
   sub(c, f, c);

   if (!divide(c, c, p))
      Error("inexact division");

   ZZ_pX cc, gg, hh, aa, bb, tt, gg1, hh1;

   conv(cc, c);
   conv(gg, g);
   conv(hh, h);
   conv(aa, a);
   conv(bb, b);

   ZZ_pXModulus GG;
   ZZ_pXModulus HH;

   build(GG, gg);
   build(HH, hh);

   ZZ_pXMultiplier AA;
   ZZ_pXMultiplier BB;

   build(AA, aa, HH);
   build(BB, bb, GG);

   rem(gg1, cc, GG);
   MulMod(gg1, gg1, BB, GG);
   
   rem(hh1, cc, HH);
   MulMod(hh1, hh1, AA, HH);

   conv(g1, gg1);
   mul(g1, g1, p);
   add(G, g, g1);

   conv(h1, hh1);
   mul(h1, h1, p);
   add(H, h, h1);

   /* lift inverses */

   ZZX t1, t2, r;

   mul(t1, a, G);
   mul(t2, b, H);
   add(t1, t1, t2);
   add(t1, t1, -1);
   negate(t1, t1);

   if (!divide(r, t1, p))
      Error("inexact division");

   ZZ_pX rr, aa1, bb1;

   conv(rr, r);
   
   rem(aa1, rr, HH);
   MulMod(aa1, aa1, AA, HH);
   rem(bb1, rr, GG);
   MulMod(bb1, bb1, BB, GG);

   ZZX a1, b1;

   conv(a1, aa1);
   mul(a1, a1, p);
   add(A, a, a1);

   conv(b1, bb1);
   mul(b1, b1, p);
   add(B, b, b1);

   Gout = G;
   Hout = H;
   Aout = A;
   Bout = B;
}

static
void HenselLift1(ZZX& Gout, ZZX& Hout, 
                const ZZX& f, const ZZX& g, const ZZX& h,
                const ZZX& a, const ZZX& b, const ZZ& p) 
{
   ZZX c, g1, h1, G, H;

   mul(c, g, h);
   sub(c, f, c);

   if (!divide(c, c, p))
      Error("inexact division");

   ZZ_pX cc, gg, hh, aa, bb, tt, gg1, hh1;

   conv(cc, c);
   conv(gg, g);
   conv(hh, h);
   conv(aa, a);
   conv(bb, b);

   ZZ_pXModulus GG;
   ZZ_pXModulus HH;

   build(GG, gg);
   build(HH, hh);

   rem(gg1, cc, GG);
   MulMod(gg1, gg1, bb, GG);
   
   rem(hh1, cc, HH);
   MulMod(hh1, hh1, aa, HH);

   conv(g1, gg1);
   mul(g1, g1, p);
   add(G, g, g1);

   conv(h1, hh1);
   mul(h1, h1, p);
   add(H, h, h1);

   Gout = G;
   Hout = H;
}

static
void BuildTree(vec_long& link, vec_ZZX& v, vec_ZZX& w,
               const vec_zz_pX& a)
{
   long k = a.length();

   if (k < 2) Error("bad arguments to BuildTree");

   vec_zz_pX V, W;

   V.SetLength(2*k-2);
   W.SetLength(2*k-2);
   link.SetLength(2*k-2);

   long i, j, s;
   long minp, mind;

   for (i = 0; i < k; i++) {
      V[i] = a[i];
      link[i] = -(i+1);
   }

   for (j = 0; j < 2*k-4; j += 2) {
      minp = j;
      mind = deg(V[j]);

      for (s = j+1; s < i; s++)
         if (deg(V[s]) < mind) {
            minp = s;
            mind = deg(V[s]);
         }

      swap(V[j], V[minp]);
      swap(link[j], link[minp]);

      minp = j+1;
      mind = deg(V[j+1]);

      for (s = j+2; s < i; s++)
         if (deg(V[s]) < mind) {
            minp = s;
            mind = deg(V[s]);
         }

      swap(V[j+1], V[minp]);
      swap(link[j+1], link[minp]);

      mul(V[i], V[j], V[j+1]);
      link[i] = j;
      i++;
   }

   zz_pX d;

   for (j = 0; j < 2*k-2; j += 2) {
      XGCD(d, W[j], W[j+1], V[j], V[j+1]);
      if (!IsOne(d))
         Error("relatively prime polynomials expected");
   }

   v.SetLength(2*k-2);
   for (j = 0; j < 2*k-2; j++)
      conv(v[j], V[j]);

   w.SetLength(2*k-2);
   for (j = 0; j < 2*k-2; j++)
      conv(w[j], W[j]);
}

static
void RecTreeLift(const vec_long& link, vec_ZZX& v, vec_ZZX& w,
                 const ZZ& p, const ZZX& f, long j, long inv)
{
   if (j < 0) return;

   if (inv)
      HenselLift(v[j], v[j+1], w[j], w[j+1],
                 f, v[j], v[j+1], w[j], w[j+1], p);
   else
      HenselLift1(v[j], v[j+1], f, v[j], v[j+1], w[j], w[j+1], p);

   RecTreeLift(link, v, w, p, v[j], link[j], inv);
   RecTreeLift(link, v, w, p, v[j+1], link[j+1], inv);
}

static
void TreeLift(const vec_long& link, vec_ZZX& v, vec_ZZX& w, 
              long e0, long e1, const ZZX& f, long inv)

// lift from p^{e0} to p^{e1}

{
   ZZ p0, p1;

   power(p0, zz_p::modulus(), e0);
   power(p1, zz_p::modulus(), e1-e0);

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

   RecTreeLift(link, v, w, p0, f, v.length()-2, inv);

   bak.restore();
} 

void MultiLift(vec_ZZX& A, const vec_zz_pX& a, const ZZX& f, long e,
               long verbose)

{
   long k = a.length();
   long i;

   if (k < 2 || e < 1) Error("MultiLift: bad args");

   if (!IsOne(LeadCoeff(f)))
      Error("MultiLift: bad args");

   for (i = 0; i < a.length(); i++)
      if (!IsOne(LeadCoeff(a[i])))
         Error("MultiLift: bad args");

   if (e == 1) {
      A.SetLength(k);
      for (i = 0; i < k; i++)
         conv(A[i], a[i]);
      return;
   }

   vec_long E;
   append(E, e);
   while (e > 1) {
      e = (e+1)/2;
      append(E, e);
   }
   long l = E.length();

   vec_ZZX v, w;
   vec_long link;

   double t;

   if (verbose) {
      cerr << "building tree...";
      t = GetTime();
   }

   BuildTree(link, v, w, a);

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


   for (i = l-1; i > 0; i--) {
      if (verbose) {
         cerr << "lifting to " << E[i-1] << "...";
         t = GetTime();
      }
      
      TreeLift(link, v, w, E[i], E[i-1], f, i != 1);

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

   A.SetLength(k);
   for (i = 0; i < 2*k-2; i++) {
      long t = link[i];
      if (t < 0)
         A[-(t+1)] = v[i];
   }
}

static
void inplace_rev(ZZX& f)
{
   long n = deg(f);
   long i, j;

   i = 0;
   j = n;
   while (i < j) {
      swap(f.rep[i], f.rep[j]);
      i++;
      j--;
   }

   f.normalize();
}

long ZZXFac_InitNumPrimes = 7;
long ZZXFac_MaxNumPrimes = 50;

static 
void RecordPattern(vec_long& pat, vec_pair_zz_pX_long& fac)
{
   long n = pat.length()-1;
   long i;

   for (i = 0; i <= n; i++)
      pat[i] = 0;

   long k = fac.length();

   for (i = 0; i < k; i++) {
      long d = fac[i].b;
      long m = deg(fac[i].a)/d;

      pat[d] = m;
   }
}

static
long NumFactors(const vec_long& pat)
{
   long n = pat.length()-1;

   long i;
   long res = 0;

   for (i = 0; i <= n; i++) 
      res += pat[i];

   return res;
}

static
void CalcPossibleDegrees(ZZ& pd, const vec_long& pat)
{
   long n = pat.length()-1;
   set(pd);

   long d, j;
   ZZ t1;

   for (d = 1; d <= n; d++) 
      for (j = 0; j < pat[d]; j++) {
         LeftShift(t1, pd, d);
         bit_or(pd, pd, t1);
      }
}

static 
void CalcPossibleDegrees(vec_ZZ& S, const vec_ZZ_pX& fac, long k)

// S[i] = possible degrees of the product of any subset of size k
//        among fac[i...], encoded as a bit vector.      

{
   long r = fac.length();

   S.SetLength(r);

   if (r == 0)
      return;

   if (k < 1 || k > r)
      Error("CalcPossibleDegrees: bad args");

   long i, l;
   ZZ old, t1;

   set(S[r-1]);
   LeftShift(S[r-1], S[r-1], deg(fac[r-1]));

   for (i = r-2; i >= 0; i--) {
      set(t1);
      LeftShift(t1, t1, deg(fac[i]));
      bit_or(S[i], t1, S[i+1]);
   }

   for (l = 2; l <= k; l++) {
      old = S[r-l];
      LeftShift(S[r-l], S[r-l+1], deg(fac[r-l]));

      for (i = r-l-1; i >= 0; i--) {
         LeftShift(t1, old, deg(fac[i]));
         old = S[i];
         bit_or(S[i], S[i+1], t1);
      }
   }
}



static
vec_zz_pX *
SmallPrimeFactorization(LocalInfoT& LocalInfo, const ZZX& f,
                            long verbose)

{
   long n = deg(f);
   long i;
   double t;

   LocalInfo.n = n;
   long& NumPrimes = LocalInfo.NumPrimes;
   NumPrimes = 0;

   LocalInfo.NumFactors = 0;

   // some sanity checking...

   if (ZZXFac_InitNumPrimes < 1 || ZZXFac_InitNumPrimes > 10000)
      Error("bad ZZXFac_InitNumPrimes");

   if (ZZXFac_MaxNumPrimes < ZZXFac_InitNumPrimes || ZZXFac_MaxNumPrimes > 10000)
      Error("bad ZZXFac_MaxNumPrimes");

   LocalInfo.p.SetLength(ZZXFac_InitNumPrimes);
   LocalInfo.pattern.SetLength(ZZXFac_InitNumPrimes);
  
   // set bits 0..n of LocalInfo.PossibleDegrees 
   SetBit(LocalInfo.PossibleDegrees, n+1);
   add(LocalInfo.PossibleDegrees, LocalInfo.PossibleDegrees, -1);

   long minr = n+1;
   long irred = 0;

   vec_pair_zz_pX_long *bestfac = 0;
   zz_pX *besth = 0;
   vec_zz_pX *spfactors = 0;
   zz_pContext bestp;
   long bestp_index;

   long maxroot = NextPowerOfTwo(deg(f))+1;

   for (; NumPrimes < ZZXFac_InitNumPrimes;) {
      long p = LocalInfo.s.next();
      if (!p) Error("out of small primes");
      if (divide(LeadCoeff(f), p)) {
         if (verbose) cerr << "skipping " << p << "\n";
         continue;
      }
      zz_p::init(p, maxroot);

      zz_pX ff, ffp, d;

      conv(ff, f);
      MakeMonic(ff);
      diff(ffp, ff);

      GCD(d, ffp, ff);
      if (!IsOne(d)) {
         if (verbose)  cerr << "skipping " << p << "\n";
         continue;
      }


      if (verbose) {
         cerr << "factoring mod " << p << "...";
         t = GetTime();
      }

      vec_pair_zz_pX_long thisfac;
      zz_pX thish; 

      SFCanZass1(thisfac, thish, ff, 0);

      LocalInfo.p[NumPrimes] = p;

      vec_long& pattern = LocalInfo.pattern[NumPrimes];
      pattern.SetLength(n+1);

      RecordPattern(pattern, thisfac);
      long r = NumFactors(pattern);
      
      if (verbose) {
         cerr << (GetTime()-t) << "\n";
         cerr << "degree sequence: ";
         for (i = 0; i <= n; i++)
            if (pattern[i]) {
               cerr << pattern[i] << "*" << i << " ";
            }
         cerr << "\n";
      }

      if (r == 1) {
         irred = 1;
         break;
      }

      // update admissibility info

      ZZ pd;

      CalcPossibleDegrees(pd, pattern);
      bit_and(LocalInfo.PossibleDegrees, LocalInfo.PossibleDegrees, pd);

      if (weight(LocalInfo.PossibleDegrees) == 2) {
         irred = 1;
         break;
      }


      if (r < minr) {
         minr = r;
         delete bestfac;
         bestfac = NTL_NEW_OP vec_pair_zz_pX_long;
         *bestfac = thisfac;
         delete besth;
         besth = NTL_NEW_OP zz_pX;
         *besth = thish;
         bestp.save();
         bestp_index = NumPrimes;
      }

      NumPrimes++;
   }