gf2exfactoring.c

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

   GF2EX y, z, t;

   z = b;
   y = a;
   clear(w);

   while (d) {
      if (d == 1) {
         if (IsZero(w)) 
            w = y;
         else {
            CompMod(w, w, z, F);
            add(w, w, y);
         }
      }
      else if ((d & 1) == 0) {
         Comp2Mod(z, t, z, y, z, F);
         add(y, t, y);
      }
      else if (IsZero(w)) {
         w = y;
         Comp2Mod(z, t, z, y, z, F);
         add(y, t, y);
      }
      else {
         Comp3Mod(z, t, w, z, y, w, z, F);
         add(w, w, y);
         add(y, t, y);
      }

      d = d >> 1;
   }
}


void PowerCompose(GF2EX& y, const GF2EX& h, long q, const GF2EXModulus& F)
{
   if (q < 0) Error("powerCompose: bad args");

   GF2EX z(INIT_SIZE, F.n);
   long sw;

   z = h;
   SetX(y);

   while (q) {
      sw = 0;

      if (q > 1) sw = 2;
      if (q & 1) {
         if (IsX(y))
            y = z;
         else
            sw = sw | 1;
      }

      switch (sw) {
      case 0:
         break;

      case 1:
         CompMod(y, y, z, F);
         break;

      case 2:
         CompMod(z, z, z, F);
         break;

      case 3:
         Comp2Mod(y, z, y, z, z, F);
         break;
      }

      q = q >> 1;
   }
}


long ProbIrredTest(const GF2EX& f, long iter)
{
   long n = deg(f);

   if (n <= 0) return 0;
   if (n == 1) return 1;

   GF2EXModulus F;

   build(F, f);

   GF2EX b, r, s;

   FrobeniusMap(b, F);

   long all_zero = 1;

   long i;

   for (i = 0; i < iter; i++) {
      random(r, n);
      TraceMap(s, r, n, F, b);

      all_zero = all_zero && IsZero(s);

      if (deg(s) > 0) return 0;
   }

   if (!all_zero || (n & 1)) return 1;

   PowerCompose(s, b, n/2, F);
   return !IsX(s);
}


long GF2EX_BlockingFactor = 10;

void DDF(vec_pair_GF2EX_long& factors, const GF2EX& ff, const GF2EX& hh, 
         long verbose)
{
   GF2EX f = ff;
   GF2EX h = hh;

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

   factors.SetLength(0);

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

   if (deg(f) == 1) {
      AddFactor(factors, f, 1, verbose);
      return;
   }

   long CompTableSize = 2*SqrRoot(deg(f)); 

   long GCDTableSize = GF2EX_BlockingFactor;

   GF2EXModulus F;
   build(F, f);

   GF2EXArgument H;

   build(H, h, F, min(CompTableSize, deg(f)));

   long i, d, limit, old_n;
   GF2EX g, X;


   vec_GF2EX tbl(INIT_SIZE, GCDTableSize);

   SetX(X);

   i = 0;
   g = h;
   d = 1;
   limit = GCDTableSize;


   while (2*d <= deg(f)) {

      old_n = deg(f);
      add(tbl[i], g, X);
      i++;
      if (i == limit) {
         ProcessTable(f, factors, F, i, tbl, d, verbose);
         i = 0;
      }

      d = d + 1;
      if (2*d <= deg(f)) {
         // we need to go further

         if (deg(f) < old_n) {
            // f has changed 

            build(F, f);
            rem(h, h, f);
            rem(g, g, f);
            build(H, h, F, min(CompTableSize, deg(f)));
         }

         CompMod(g, g, H, F);
      }
   }

   ProcessTable(f, factors, F, i, tbl, d-1, verbose);

   if (!IsOne(f)) AddFactor(factors, f, deg(f), verbose);
}



void RootEDF(vec_GF2EX& factors, const GF2EX& f, long verbose)
{
   vec_GF2E roots;
   double t;

   if (verbose) { cerr << "finding roots..."; t = GetTime(); }
   FindRoots(roots, f);
   if (verbose) { cerr << (GetTime()-t) << "\n"; }

   long r = roots.length();
   factors.SetLength(r);
   for (long j = 0; j < r; j++) {
      SetX(factors[j]);
      add(factors[j], factors[j], roots[j]);
   }
}

static
void EDFSplit(vec_GF2EX& v, const GF2EX& f, const GF2EX& b, long d)
{
   GF2EX a, g, h;
   GF2EXModulus F;
   vec_GF2E roots;
   
   build(F, f);
   long n = F.n;
   long r = n/d;
   random(a, n);
   TraceMap(g, a, d, F, b);
   MinPolyMod(h, g, F, r);
   FindRoots(roots, h);
   FindFactors(v, f, g, roots);
}

static
void RecEDF(vec_GF2EX& factors, const GF2EX& f, const GF2EX& b, long d,
            long verbose)
{
   vec_GF2EX v;
   long i;
   GF2EX bb;

   if (verbose) cerr << "+";

   EDFSplit(v, f, b, d);
   for (i = 0; i < v.length(); i++) {
      if (deg(v[i]) == d) {
         append(factors, v[i]);
      }
      else {
         GF2EX bb;
         rem(bb, b, v[i]);
         RecEDF(factors, v[i], bb, d, verbose);
      }
   }
}
         

void EDF(vec_GF2EX& factors, const GF2EX& ff, const GF2EX& bb,
         long d, long verbose)

{
   GF2EX f = ff;
   GF2EX b = bb;

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

   long n = deg(f);
   long r = n/d;

   if (r == 0) {
      factors.SetLength(0);
      return;
   }

   if (r == 1) {
      factors.SetLength(1);
      factors[0] = f;
      return;
   }

   if (d == 1) {
      RootEDF(factors, f, verbose);
      return;
   }

   
   double t;
   if (verbose) { 
      cerr << "computing EDF(" << d << "," << r << ")..."; 
      t = GetTime(); 
   }

   factors.SetLength(0);

   RecEDF(factors, f, b, d, verbose);

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


void SFCanZass(vec_GF2EX& factors, const GF2EX& ff, long verbose)
{
   GF2EX f = ff;

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

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

   if (deg(f) == 1) {
      factors.SetLength(1);
      factors[0] = f;
      return;
   }

   factors.SetLength(0);

   double t;

   
   GF2EXModulus F;
   build(F, f);

   GF2EX h;

   if (verbose) { cerr << "computing X^p..."; t = GetTime(); }
   FrobeniusMap(h, F);
   if (verbose) { cerr << (GetTime()-t) << "\n"; }

   vec_pair_GF2EX_long u;
   if (verbose) { cerr << "computing DDF..."; t = GetTime(); }
   NewDDF(u, f, h, verbose);
   if (verbose) { 
      t = GetTime()-t; 
      cerr << "DDF time: " << t << "\n";
   }

   GF2EX hh;
   vec_GF2EX v;

   long i;
   for (i = 0; i < u.length(); i++) {
      const GF2EX& g = u[i].a;
      long d = u[i].b;
      long r = deg(g)/d;

      if (r == 1) {
         // g is already irreducible

         append(factors, g);
      }
      else {
         // must perform EDF

         if (d == 1) {
            // root finding
            RootEDF(v, g, verbose);
            append(factors, v);
         }
         else {
            // general case
            rem(hh, h, g);
            EDF(v, g, hh, d, verbose);
            append(factors, v);
         }
      }
   }
}
   
void CanZass(vec_pair_GF2EX_long& factors, const GF2EX& f, long verbose)
{
   if (!IsOne(LeadCoeff(f)))
      Error("CanZass: bad args");

   double t;
   vec_pair_GF2EX_long sfd;
   vec_GF2EX x;

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

   factors.SetLength(0);

   long i, j;

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

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

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

void mul(GF2EX& f, const vec_pair_GF2EX_long& v)
{
   long i, j, n;

   n = 0;
   for (i = 0; i < v.length(); i++)
      n += v[i].b*deg(v[i].a);

   GF2EX g(INIT_SIZE, n+1);

   set(g);
   for (i = 0; i < v.length(); i++)
      for (j = 0; j < v[i].b; j++) {
         mul(g, g, v[i].a);
      }

   f = g;
}




static
long BaseCase(const GF2EX& h, long q, long a, const GF2EXModulus& F)
{
   long b, e;
   GF2EX lh(INIT_SIZE, F.n);

   lh = h;
   b = 1;
   e = 0;
   while (e < a-1 && !IsX(lh)) {
      e++;
      b *= q;
      PowerCompose(lh, lh, q, F);
   }

   if (!IsX(lh)) b *= q;

   return b;
}



static
void TandemPowerCompose(GF2EX& y1, GF2EX& y2, const GF2EX& h, 
                        long q1, long q2, const GF2EXModulus& F)
{
   GF2EX z(INIT_SIZE, F.n);
   long sw;

   z = h;
   SetX(y1);
   SetX(y2);

   while (q1 || q2) {
      sw = 0;

      if (q1 > 1 || q2 > 1) sw = 4;

      if (q1 & 1) {
         if (IsX(y1))
            y1 = z;
         else
            sw = sw | 2;
      }

      if (q2 & 1) {
         if (IsX(y2))
            y2 = z;
         else
            sw = sw | 1;
      }

      switch (sw) {
      case 0:
         break;

      case 1:
         CompMod(y2, y2, z, F);
         break;

      case 2:
         CompMod(y1, y1, z, F);
         break;

      case 3:
         Comp2Mod(y1, y2, y1, y2, z, F);
         break;

      case 4:
         CompMod(z, z, z, F);
         break;

      case 5:
         Comp2Mod(z, y2, z, y2, z, F);
         break;

      case 6:
         Comp2Mod(z, y1, z, y1, z, F);
         break;

      case 7:
         Comp3Mod(z, y1, y2, z, y1, y2, z, F);
         break;
      }

      q1 = q1 >> 1;
      q2 = q2 >> 1;
   }
}



static
long RecComputeDegree(long u, const GF2EX& h, const GF2EXModulus& F,
                      FacVec& fvec)
{
   if (IsX(h)) return 1;

   if (fvec[u].link == -1) return BaseCase(h, fvec[u].q, fvec[u].a, F);

   GF2EX h1, h2;
   long q1, q2, r1, r2;

   q1 = fvec[fvec[u].link].val; 
   q2 = fvec[fvec[u].link+1].val;

   TandemPowerCompose(h1, h2, h, q1, q2, F);
   r1 = RecComputeDegree(fvec[u].link, h2, F, fvec);
   r2 = RecComputeDegree(fvec[u].link+1, h1, F, fvec);
   return r1*r2;
}

   


long RecComputeDegree(const GF2EX& h, const GF2EXModulus& F)
   // f = F.f is assumed to be an "equal degree" polynomial
   // h = X^p mod f
   // the common degree of the irreducible factors of f is computed
{
   if (F.n == 1 || IsX(h)) 
      return 1;

   FacVec fvec;

   FactorInt(fvec, F.n);

   return RecComputeDegree(fvec.length()-1, h, F, fvec);
}


void FindRoot(GF2E& root, const GF2EX& ff)
// finds a root of ff.
// assumes that ff is monic and splits into distinct linear factors

{
   GF2EXModulus F;
   GF2EX h, h1, f;
   GF2E r;

   f = ff;
   
   if (!IsOne(LeadCoeff(f)))
      Error("FindRoot: bad args");

   if (deg(f) == 0)
      Error("FindRoot: bad args");


   while (deg(f) > 1) {
      build(F, f);
      random(r);
      clear(h);
      SetCoeff(h, 1, r);
      TraceMap(h, h, F);
      GCD(h, h, f);
      if (deg(h) > 0 && deg(h) < deg(f)) {
         if (deg(h) > deg(f)/2)
            div(f, f, h);
         else
            f = h;
      }
   }
 
   root = ConstTerm(f);
}


static
long power(long a, long e)
{
   long i, res;

   res = 1;
   for (i = 1; i <= e; i++)
      res = res * a;

   return res;
}


static
long IrredBaseCase(const GF2EX& h, long q, long a, const GF2EXModulus& F)
{
   long e;
   GF2EX X, s, d;

   e = power(q, a-1);
   PowerCompose(s, h, e, F);
   SetX(X);
   add(s, s, X);
   GCD(d, F.f, s);
   return IsOne(d);
}


static
long RecIrredTest(long u, const GF2EX& h, const GF2EXModulus& F,
                 const FacVec& fvec)
{
   long  q1, q2;
   GF2EX h1, h2;

   if (IsX(h)) return 0;

   if (fvec[u].link == -1) {
      return IrredBaseCase(h, fvec[u].q, fvec[u].a, F);
   }


   q1 = fvec[fvec[u].link].val; 
   q2 = fvec[fvec[u].link+1].val;

   TandemPowerCompose(h1, h2, h, q1, q2, F);
   return RecIrredTest(fvec[u].link, h2, F, fvec) 
          && RecIrredTest(fvec[u].link+1, h1, F, fvec);
}

long DetIrredTest(const GF2EX& f)
{
   if (deg(f) <= 0) return 0;
   if (deg(f) == 1) return 1;

   GF2EXModulus F;

   build(F, f);
   
   GF2EX h;

   FrobeniusMap(h, F);

   GF2EX s;
   PowerCompose(s, h, F.n, F);
   if (!IsX(s)) return 0;

   FacVec fvec;

   FactorInt(fvec, F.n);

   return RecIrredTest(fvec.length()-1, h, F, fvec);
}



long IterIrredTest(const GF2EX& f)
{
   if (deg(f) <= 0) return 0;
   if (deg(f) == 1) return 1;

   GF2EXModulus F;

   build(F, f);
   
   GF2EX h;

   FrobeniusMap(h, F);

   long CompTableSize = 2*SqrRoot(deg(f));

   GF2EXArgument H;

   build(H, h, F, CompTableSize);

   long i, d, limit, limit_sqr;
   GF2EX g, X, t, prod;


   SetX(X);

   i = 0;
   g = h;
   d = 1;
   limit = 2;
   limit_sqr = limit*limit;

   set(prod);


   while (2*d <= deg(f)) {
      add(t, g, X);
      MulMod(prod, prod, t, F);
      i++;
      if (i == limit_sqr) {
         GCD(t, f, prod);
         if (!IsOne(t)) return 0;

         set(prod);
         limit++;
         limit_sqr = limit*limit;
         i = 0;
      }

      d = d + 1;
      if (2*d <= deg(f)) {
         CompMod(g, g, H, F);
      }
   }

   if (i > 0) {
      GCD(t, f, prod);
      if (!IsOne(t)) return 0;
   }

   return 1;
}

static
void MulByXPlusY(vec_GF2EX& h, const GF2EX& f, const GF2EX& g)
// h represents the bivariate polynomial h[0] + h[1]*Y + ... + h[n-1]*Y^k,
// where the h[i]'s are polynomials in X, each of degree < deg(f),
// and k < deg(g).