gf2exfactoring.c

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

// h is replaced by the bivariate polynomial h*(X+Y) (mod f(X), g(Y)).

{
   long n = deg(g);
   long k = h.length()-1;

   if (k < 0) return;

   if (k < n-1) {
      h.SetLength(k+2);
      h[k+1] = h[k];
      for (long i = k; i >= 1; i--) {
         MulByXMod(h[i], h[i], f);
         add(h[i], h[i], h[i-1]);
      }
      MulByXMod(h[0], h[0], f);
   }
   else {
      GF2EX b, t;

      b = h[n-1];
      for (long i = n-1; i >= 1; i--) {
         mul(t, b, g.rep[i]);
         MulByXMod(h[i], h[i], f);
         add(h[i], h[i], h[i-1]);
         add(h[i], h[i], t);
      }
      mul(t, b, g.rep[0]);
      MulByXMod(h[0], h[0], f);
      add(h[0], h[0], t);
   }

   // normalize

   k = h.length()-1;
   while (k >= 0 && IsZero(h[k])) k--;
   h.SetLength(k+1);
}


static
void IrredCombine(GF2EX& x, const GF2EX& f, const GF2EX& g)
{
   if (deg(f) < deg(g)) {
      IrredCombine(x, g, f);
      return;
   }

   // deg(f) >= deg(g)...not necessary, but maybe a little more
   //                    time & space efficient

   long df = deg(f);
   long dg = deg(g);
   long m = df*dg;

   vec_GF2EX h(INIT_SIZE, dg);

   long i;
   for (i = 0; i < dg; i++) h[i].SetMaxLength(df);

   h.SetLength(1);
   set(h[0]);

   vec_GF2E a;

   a.SetLength(2*m);

   for (i = 0; i < 2*m; i++) {
      a[i] = ConstTerm(h[0]);
      if (i < 2*m-1)
         MulByXPlusY(h, f, g);
   }

   MinPolySeq(x, a, m);
}


static
void BuildPrimePowerIrred(GF2EX& f, long q, long e)
{
   long n = power(q, e);

   do {
      random(f, n);
      SetCoeff(f, n);
   } while (!IterIrredTest(f));
}

static
void RecBuildIrred(GF2EX& f, long u, const FacVec& fvec)
{
   if (fvec[u].link == -1)
      BuildPrimePowerIrred(f, fvec[u].q, fvec[u].a);
   else {
      GF2EX g, h;
      RecBuildIrred(g, fvec[u].link, fvec);
      RecBuildIrred(h, fvec[u].link+1, fvec);
      IrredCombine(f, g, h);
   }
}


void BuildIrred(GF2EX& f, long n)
{
   if (n <= 0)
      Error("BuildIrred: n must be positive");

   if (n >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in BuildIrred");

   if (n == 1) {
      SetX(f);
      return;
   }

   FacVec fvec;

   FactorInt(fvec, n);

   RecBuildIrred(f, fvec.length()-1, fvec);
}



#if 0
void BuildIrred(GF2EX& f, long n)
{
   if (n <= 0)
      Error("BuildIrred: n must be positive");

   if (n >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in BuildIrred");

   if (n == 1) {
      SetX(f);
      return;
   }

   GF2EX g;

   do {
      random(g, n);
      SetCoeff(g, n);
   } while (!IterIrredTest(g));

   f = g;

}
#endif



void BuildRandomIrred(GF2EX& f, const GF2EX& g)
{
   GF2EXModulus G;
   GF2EX h, ff;

   build(G, g);
   do {
      random(h, deg(g));
      IrredPolyMod(ff, h, G);
   } while (deg(ff) < deg(g));

   f = ff;
}


/************* NEW DDF ****************/

long GF2EX_GCDTableSize = 4;
char GF2EX_stem[256] = "";

double GF2EXFileThresh = 256;

static vec_GF2EX BabyStepFile;
static vec_GF2EX GiantStepFile;

static long use_files;


static
double CalcTableSize(long n, long k)
{
   double sz = GF2E::WordLength()+4;
   sz = sz * sizeof(_ntl_ulong);
   sz = sz * n;
   sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_GF2E);
   sz = sz * k;
   sz = sz/1024;
   return sz;
}



static
void GenerateBabySteps(GF2EX& h1, const GF2EX& f, const GF2EX& h, long k,
                       long verbose)

{
   double t;

   if (verbose) { cerr << "generating baby steps..."; t = GetTime(); }

   GF2EXModulus F;
   build(F, f);

   GF2EXArgument H;

#if 0
   double n2 = sqrt(double(F.n));
   double n4 = sqrt(n2);
   double n34 = n2*n4;
   long sz = long(ceil(n34/sqrt(sqrt(2.0))));
#else
   long sz = 2*SqrRoot(F.n);
#endif

   build(H, h, F, sz);


   h1 = h;

   long i;

   long HexOutput = GF2X::HexOutput;
   GF2X::HexOutput = 1;

   if (!use_files) {
      BabyStepFile.kill();
      BabyStepFile.SetLength(k-1);
   }

   for (i = 1; i <= k-1; i++) {
      if (use_files) {
         ofstream s;
         OpenWrite(s, FileName(GF2EX_stem, "baby", i));
         s << h1 << "\n";
         s.close();
      }
      else
         BabyStepFile(i) = h1;

      CompMod(h1, h1, H, F);
      if (verbose) cerr << "+";
   }

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

   GF2X::HexOutput = HexOutput;
}


static
void GenerateGiantSteps(const GF2EX& f, const GF2EX& h, long l, long verbose)
{

   double t;

   if (verbose) { cerr << "generating giant steps..."; t = GetTime(); }

   GF2EXModulus F;
   build(F, f);

   GF2EXArgument H;

#if 0
   double n2 = sqrt(double(F.n));
   double n4 = sqrt(n2);
   double n34 = n2*n4;
   long sz = long(ceil(n34/sqrt(sqrt(2.0))));
#else
   long sz = 2*SqrRoot(F.n);
#endif

   build(H, h, F, sz);

   GF2EX h1;

   h1 = h;

   long i;

   long HexOutput = GF2X::HexOutput; 
   GF2X::HexOutput = 1;

   if (!use_files) {
      GiantStepFile.kill();
      GiantStepFile.SetLength(l);
   }

   for (i = 1; i <= l-1; i++) {
      if (use_files) {
         ofstream s;
         OpenWrite(s, FileName(GF2EX_stem, "giant", i));
         s << h1 << "\n";
         s.close();
      }
      else
         GiantStepFile(i) = h1;

      CompMod(h1, h1, H, F);
      if (verbose) cerr << "+";
   }

   if (use_files) {
      ofstream s;
      OpenWrite(s, FileName(GF2EX_stem, "giant", i));
      s << h1 << "\n";
      s.close();
   }
   else
      GiantStepFile(i) = h1;

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

   GF2X::HexOutput = HexOutput;
}

static
void FileCleanup(long k, long l)
{
   if (use_files) {
      long i;
   
      for (i = 1; i <= k-1; i++)
         remove(FileName(GF2EX_stem, "baby", i));
   
      for (i = 1; i <= l; i++)
         remove(FileName(GF2EX_stem, "giant", i));
   }
   else {
      BabyStepFile.kill();
      GiantStepFile.kill();
   }
}


static
void NewAddFactor(vec_pair_GF2EX_long& u, const GF2EX& g, long m, long verbose)
{
   long len = u.length();

   u.SetLength(len+1);
   u[len].a = g;
   u[len].b = m;

   if (verbose) {
      cerr << "split " << m << " " << deg(g) << "\n";
   }
}

   


static
void NewProcessTable(vec_pair_GF2EX_long& u, GF2EX& f, const GF2EXModulus& F,
                     vec_GF2EX& buf, long size, long StartInterval,
                     long IntervalLength, long verbose)

{
   if (size == 0) return;

   GF2EX& g = buf[size-1];

   long i;

   for (i = 0; i < size-1; i++)
      MulMod(g, g, buf[i], F);

   GCD(g, f, g);

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

   div(f, f, g);

   long d = (StartInterval-1)*IntervalLength + 1;
   i = 0;
   long interval = StartInterval;

   while (i < size-1 && 2*d <= deg(g)) {
      GCD(buf[i], buf[i], g);
      if (deg(buf[i]) > 0) {
         NewAddFactor(u, buf[i], interval, verbose);
         div(g, g, buf[i]);
      }

      i++;
      interval++;
      d += IntervalLength;
   }

   if (deg(g) > 0) {
      if (i == size-1)
         NewAddFactor(u, g, interval, verbose);
      else
         NewAddFactor(u, g, (deg(g)+IntervalLength-1)/IntervalLength, verbose);
   }
}


static
void FetchGiantStep(GF2EX& g, long gs, const GF2EXModulus& F)
{
   if (use_files) {
      ifstream s;
   
      OpenRead(s, FileName(GF2EX_stem, "giant", gs));
   
      s >> g;
      s.close();
   }
   else
      g = GiantStepFile(gs);

   rem(g, g, F);
}


static
void FetchBabySteps(vec_GF2EX& v, long k)
{
   v.SetLength(k);

   SetX(v[0]);

   long i;
   for (i = 1; i <= k-1; i++) {
      if (use_files) {
         ifstream s;
         OpenRead(s, FileName(GF2EX_stem, "baby", i));
         s >> v[i];
         s.close();
      }
      else
         v[i] = BabyStepFile(i);
   }
}
      


static
void GiantRefine(vec_pair_GF2EX_long& u, const GF2EX& ff, long k, long l,
                 long verbose)

{
   double t;

   if (verbose) {
      cerr << "giant refine...";
      t = GetTime();
   }

   u.SetLength(0);

   vec_GF2EX BabyStep;

   FetchBabySteps(BabyStep, k);

   vec_GF2EX buf(INIT_SIZE, GF2EX_GCDTableSize);

   GF2EX f;
   f = ff;

   GF2EXModulus F;
   build(F, f);

   GF2EX g;
   GF2EX h;

   long size = 0;

   long first_gs;

   long d = 1;

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

      long old_n = deg(f);

      long gs = (d+k-1)/k;
      long bs = gs*k - d;

      if (bs == k-1) {
         size++;
         if (size == 1) first_gs = gs;
         FetchGiantStep(g, gs, F);
         add(buf[size-1], g, BabyStep[bs]);
      }
      else {
         add(h, g, BabyStep[bs]);
         MulMod(buf[size-1], buf[size-1], h, F);
      }

      if (verbose && bs == 0) cerr << "+";

      if (size == GF2EX_GCDTableSize && bs == 0) {
         NewProcessTable(u, f, F, buf, size, first_gs, k, verbose);
         if (verbose) cerr << "*";
         size = 0;
      }

      d++;

      if (2*d <= deg(f) && deg(f) < old_n) {
         build(F, f);

         long i;
         for (i = 1; i <= k-1; i++) 
            rem(BabyStep[i], BabyStep[i], F);
      }
   }

   if (size > 0) {
      NewProcessTable(u, f, F, buf, size, first_gs, k, verbose);
      if (verbose) cerr << "*";
   }

   if (deg(f) > 0) 
      NewAddFactor(u, f, 0, verbose);

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


static
void IntervalRefine(vec_pair_GF2EX_long& factors, const GF2EX& ff,
                    long k, long gs, const vec_GF2EX& BabyStep, long verbose)

{
   vec_GF2EX buf(INIT_SIZE, GF2EX_GCDTableSize);

   GF2EX f;
   f = ff;

   GF2EXModulus F;
   build(F, f);

   GF2EX g;

   FetchGiantStep(g, gs, F);

   long size = 0;

   long first_d;

   long d = (gs-1)*k + 1;
   long bs = k-1;

   while (bs >= 0 && 2*d <= deg(f)) {

      long old_n = deg(f);

      if (size == 0) first_d = d;
      rem(buf[size], BabyStep[bs], F);
      add(buf[size], buf[size], g);
      size++;

      if (size == GF2EX_GCDTableSize) {
         NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose);
         size = 0;
      }

      d++;
      bs--;

      if (bs >= 0 && 2*d <= deg(f) && deg(f) < old_n) {
         build(F, f);
         rem(g, g, F);
      }
   }

   NewProcessTable(factors, f, F, buf, size, first_d, 1, verbose);

   if (deg(f) > 0) 
      NewAddFactor(factors, f, deg(f), verbose);
}
   



static
void BabyRefine(vec_pair_GF2EX_long& factors, const vec_pair_GF2EX_long& u,
                long k, long l, long verbose)

{
   double t;

   if (verbose) {
      cerr << "baby refine...";
      t = GetTime();
   }

   factors.SetLength(0);

   vec_GF2EX BabyStep;

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

      if (gs == 0 || 2*((gs-1)*k+1) > deg(g))
         NewAddFactor(factors, g, deg(g), verbose);
      else {
         if (BabyStep.length() == 0)
            FetchBabySteps(BabyStep, k);
         IntervalRefine(factors, g, k, gs, BabyStep, verbose);
      }
   }

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

      
      

      

void NewDDF(vec_pair_GF2EX_long& factors,
            const GF2EX& f,
            const GF2EX& h,
            long verbose)

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

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

   if (deg(f) == 1) {
      factors.SetLength(0);
      append(factors, cons(f, 1));
      return;
   }

   if (!GF2EX_stem[0])
      sprintf(GF2EX_stem, "ddf-%ld", RandomBnd(10000));
      
   long B = deg(f)/2;
   long k = SqrRoot(B);
   long l = (B+k-1)/k;

   GF2EX h1;

   if (CalcTableSize(deg(f), k + l - 1) > GF2EXFileThresh)
      use_files = 1;
   else
      use_files = 0;

   GenerateBabySteps(h1, f, h, k, verbose);

   GenerateGiantSteps(f, h1, l, verbose);

   vec_pair_GF2EX_long u;
   GiantRefine(u, f, k, l, verbose);
   BabyRefine(factors, u, k, l, verbose);

   FileCleanup(k, l);
}

long IterComputeDegree(const GF2EX& h, const GF2EXModulus& F)
{
   long n = deg(F);

   if (n == 1 || IsX(h)) return 1;

   long B = n/2;
   long k = SqrRoot(B);
   long l = (B+k-1)/k;


   GF2EXArgument H;

#if 0
   double n2 = sqrt(double(n));
   double n4 = sqrt(n2);
   double n34 = n2*n4;
   long sz = long(ceil(n34/sqrt(sqrt(2.0))));
#else
   long sz = 2*SqrRoot(F.n);
#endif

   build(H, h, F, sz);

   GF2EX h1;
   h1 = h;

   vec_GF2EX baby;
   baby.SetLength(k);

   SetX(baby[0]);

   long i;

   for (i = 1; i <= k-1; i++) {
      baby[i] = h1;
      CompMod(h1, h1, H, F);
      if (IsX(h1)) return i+1;
   }

   build(H, h1, F, sz);

   long j;

   for (j = 2; j <= l; j++) {
      CompMod(h1, h1, H, F);

      for (i = k-1; i >= 0; i--) {
         if (h1 == baby[i])
            return j*k-i;
      }
   }

   return n;
}

NTL_END_IMPL