zz_px1.cpp

数值算法库for Windows

         add(a[i], t, a[i-1]);
      }
      mul(a[0], a[0], b);
   }
} 

void mul(ZZ_p* x, const ZZ_p* a, const ZZ_p* b, long n)
{
   static ZZ t, accum;

   long i, j, jmin, jmax;

   long d = 2*n-1;

   for (i = 0; i <= d; i++) {
      jmin = max(0, i-(n-1));
      jmax = min(n-1, i);
      clear(accum);
      for (j = jmin; j <= jmax; j++) {
         mul(t, rep(a[j]), rep(b[i-j]));
         add(accum, accum, t);
      }
      if (i >= n) {
         add(accum, accum, rep(a[i-n]));
         add(accum, accum, rep(b[i-n]));
      }

      conv(x[i], accum);
   }
}


void BuildFromRoots(ZZ_pX& x, const vec_ZZ_p& a)
{
   long n = a.length();

   if (n == 0) {
      set(x);
      return;
   }

   long k0 = NextPowerOfTwo(NTL_ZZ_pX_FFT_CROSSOVER);
   long crossover = 1L << k0;

   if (n <= crossover) {
      x.rep.SetMaxLength(n+1);
      x.rep = a;
      IterBuild(&x.rep[0], n);
      x.rep.SetLength(n+1);
      SetCoeff(x, n);
      return;
   }

   long k = NextPowerOfTwo(n);

   long m = 1L << k;
   long i, j;
   long l, width;

   ZZ_pX b(INIT_SIZE, m+1);

   b.rep = a;
   b.rep.SetLength(m+1);
   for (i = n; i < m; i++)
      clear(b.rep[i]);

   set(b.rep[m]);
   
   FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);


   ZZ_p t1, one;
   set(one);

   vec_ZZ_p G(INIT_SIZE, crossover), H(INIT_SIZE, crossover);
   ZZ_p *g = G.elts();
   ZZ_p *h = H.elts();
   ZZ_p *tmp;
   
   for (i = 0; i < m; i+= crossover) {
      for (j = 0; j < crossover; j++)
         negate(g[j], b.rep[i+j]);

      if (k0 > 0) {
         for (j = 0; j < crossover; j+=2) {
            mul(t1, g[j], g[j+1]);
            add(g[j+1], g[j], g[j+1]);
            g[j] = t1;
         }
      }
   
      for (l = 1; l < k0; l++) {
         width = 1L << l;

         for (j = 0; j < crossover; j += 2*width)
            mul(&h[j], &g[j], &g[j+width], width);
      
         tmp = g; g = h; h = tmp;
      }

      for (j = 0; j < crossover; j++)
         b.rep[i+j] = g[j];
   }

   for (l = k0; l < k; l++) {
      width = 1L << l;
      for (i = 0; i < m; i += 2*width) {
         t1 = b.rep[i+width];
         set(b.rep[i+width]);
         ToFFTRep(R1, b, l+1, i, i+width);
         b.rep[i+width] = t1;
         t1 = b.rep[i+2*width];
         set(b.rep[i+2*width]);
         ToFFTRep(R2, b, l+1, i+width, i+2*width);
         b.rep[i+2*width] = t1;
         mul(R1, R1, R2);
         FromFFTRep(&b.rep[i], R1, 0, 2*width-1);
         sub(b.rep[i], b.rep[i], one);
      }
   }

   x.rep.SetLength(n+1);
   long delta = m-n;
   for (i = 0; i <= n; i++)
     x.rep[i] = b.rep[i+delta];

   // no need to normalize
}



void eval(ZZ_p& b, const ZZ_pX& f, const ZZ_p& a)
// does a Horner evaluation
{
   ZZ_p acc;
   long i;

   clear(acc);
   for (i = deg(f); i >= 0; i--) {
      mul(acc, acc, a);
      add(acc, acc, f.rep[i]);
   }

   b = acc;
}



void eval(vec_ZZ_p& b, const ZZ_pX& f, const vec_ZZ_p& a)
// naive algorithm:  repeats Horner
{
   if (&b == &f.rep) {
      vec_ZZ_p bb;
      eval(bb, f, a);
      b = bb;
      return;
   }

   long m = a.length();
   b.SetLength(m);
   long i;
   for (i = 0; i < m; i++) 
      eval(b[i], f, a[i]);
}




void interpolate(ZZ_pX& f, const vec_ZZ_p& a, const vec_ZZ_p& b)
{
   long m = a.length();
   if (b.length() != m) Error("interpolate: vector length mismatch");

   if (m == 0) {
      clear(f);
      return;
   }

   vec_ZZ_p prod;
   prod = a;

   ZZ_p t1, t2;

   long k, i;

   vec_ZZ_p res;
   res.SetLength(m);

   for (k = 0; k < m; k++) {

      const ZZ_p& aa = a[k];

      set(t1);
      for (i = k-1; i >= 0; i--) {
         mul(t1, t1, aa);
         add(t1, t1, prod[i]);
      }

      clear(t2);
      for (i = k-1; i >= 0; i--) {
         mul(t2, t2, aa);
         add(t2, t2, res[i]);
      }


      inv(t1, t1);
      sub(t2, b[k], t2);
      mul(t1, t1, t2);

      for (i = 0; i < k; i++) {
         mul(t2, prod[i], t1);
         add(res[i], res[i], t2);
      }

      res[k] = t1;

      if (k < m-1) {
         if (k == 0)
            negate(prod[0], prod[0]);
         else {
            negate(t1, a[k]);
            add(prod[k], t1, prod[k-1]);
            for (i = k-1; i >= 1; i--) {
               mul(t2, prod[i], t1);
               add(prod[i], t2, prod[i-1]);
            }
            mul(prod[0], prod[0], t1);
         }
      }
   }

   while (m > 0 && IsZero(res[m-1])) m--; 
   res.SetLength(m);
   f.rep = res;
}

NTL_vector_impl(ZZ_pX,vec_ZZ_pX)

NTL_eq_vector_impl(ZZ_pX,vec_ZZ_pX)

NTL_io_vector_impl(ZZ_pX,vec_ZZ_pX)



   
void InnerProduct(ZZ_pX& x, const vec_ZZ_p& v, long low, long high, 
                   const vec_ZZ_pX& H, long n, ZZVec& t)
{
   static ZZ s;
   long i, j;

   for (j = 0; j < n; j++)
      clear(t[j]);

   high = min(high, v.length()-1);
   for (i = low; i <= high; i++) {
      const vec_ZZ_p& h = H[i-low].rep;
      long m = h.length();
      const ZZ& w = rep(v[i]);

      for (j = 0; j < m; j++) {
         mul(s, w, rep(h[j]));
         add(t[j], t[j], s);
      }
   }

   x.rep.SetLength(n);
   for (j = 0; j < n; j++)
      conv(x.rep[j], t[j]);
   x.normalize();
}


void CompMod(ZZ_pX& x, const ZZ_pX& g, const ZZ_pXArgument& A, 
             const ZZ_pXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   ZZ_pX s, t;
   ZZVec scratch(F.n, ZZ_pInfo->ExtendedModulusSize);

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   ZZ_pXMultiplier M;
   build(M, A.H[m], F);

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }

   x = t;
}


void build(ZZ_pXArgument& A, const ZZ_pX& h, const ZZ_pXModulus& F, long m)
{
   if (m <= 0 || deg(h) >= F.n) Error("build: bad args");

   if (m > F.n) m = F.n;

   long i;

   if (ZZ_pXArgBound > 0) {
      double sz = ZZ_p::storage();
      sz = sz*F.n;
      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_ZZ_p);
      sz = sz/1024;
      m = min(m, long(ZZ_pXArgBound/sz));
      m = max(m, 1);
   }

   ZZ_pXMultiplier M;

   build(M, h, F);

   A.H.SetLength(m+1);

   set(A.H[0]);
   A.H[1] = h;
   for (i = 2; i <= m; i++) 
      MulMod(A.H[i], A.H[i-1], M, F);
}




long ZZ_pXArgBound = 0;


void CompMod(ZZ_pX& x, const ZZ_pX& g, const ZZ_pX& h, const ZZ_pXModulus& F)
   // x = g(h) mod f
{
   long m = SqrRoot(g.rep.length());

   if (m == 0) {
      clear(x);
      return;
   }

   ZZ_pXArgument A;

   build(A, h, F, m);

   CompMod(x, g, A, F);
}




void Comp2Mod(ZZ_pX& x1, ZZ_pX& x2, const ZZ_pX& g1, const ZZ_pX& g2,
              const ZZ_pX& h, const ZZ_pXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length());

   if (m == 0) {
      clear(x1);
      clear(x2);
      return;
   }

   ZZ_pXArgument A;

   build(A, h, F, m);

   ZZ_pX xx1, xx2;

   CompMod(xx1, g1, A, F);
   CompMod(xx2, g2, A, F);

   x1 = xx1;
   x2 = xx2;
}

void Comp3Mod(ZZ_pX& x1, ZZ_pX& x2, ZZ_pX& x3, 
              const ZZ_pX& g1, const ZZ_pX& g2, const ZZ_pX& g3,
              const ZZ_pX& h, const ZZ_pXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());

   if (m == 0) {
      clear(x1);
      clear(x2);
      clear(x3);
      return;
   }

   ZZ_pXArgument A;

   build(A, h, F, m);

   ZZ_pX xx1, xx2, xx3;

   CompMod(xx1, g1, A, F);
   CompMod(xx2, g2, A, F);
   CompMod(xx3, g3, A, F);

   x1 = xx1;
   x2 = xx2;
   x3 = xx3;
}


static void StripZeroes(vec_ZZ_p& x)
{
   long n = x.length();
   while (n > 0 && IsZero(x[n-1]))
      n--;
   x.SetLength(n);
}


void PlainUpdateMap(vec_ZZ_p& xx, const vec_ZZ_p& a, 
                    const ZZ_pX& b, const ZZ_pX& f)
{
   long n = deg(f);
   long i, m;

   if (IsZero(b)) {
      xx.SetLength(0);
      return;
   }

   m = n-1 - deg(b);

   vec_ZZ_p x(INIT_SIZE, n);

   for (i = 0; i <= m; i++)
      InnerProduct(x[i], a, b.rep, i);

   if (deg(b) != 0) {
      ZZ_pX c(INIT_SIZE, n);
      LeftShift(c, b, m);

      for (i = m+1; i < n; i++) {
         MulByXMod(c, c, f);
         InnerProduct(x[i], a, c.rep);
      }
   }

   xx = x;
}
   

void UpdateMap(vec_ZZ_p& x, const vec_ZZ_p& aa, 
               const ZZ_pXMultiplier& B, const ZZ_pXModulus& F)
{
   long n = F.n;
   long i;


   vec_ZZ_p a;
   a = aa;
   StripZeroes(a);

   if (a.length() > n) Error("UpdateMap: bad args");

   if (!B.UseFFT) {
      PlainUpdateMap(x, a, B.b, F.f);
      StripZeroes(x);
      return;
   }

   FFTRep R1(INIT_SIZE, F.k), R2(INIT_SIZE, F.l);
   vec_ZZ_p V1(INIT_SIZE, n);


   RevToFFTRep(R1, a, F.k, 0, a.length()-1, 0);
   mul(R2, R1, F.FRep);
   RevFromFFTRep(V1, R2, 0, n-2);
   for (i = 0; i <= n-2; i++)  negate(V1[i], V1[i]);
   RevToFFTRep(R2, V1, F.l, 0, n-2, n-1);
   mul(R2, R2, B.B1);
   mul(R1, R1, B.B2);

   AddExpand(R2, R1);
   RevFromFFTRep(x, R2, 0, n-1);
   StripZeroes(x);
}

   

void ProjectPowers(vec_ZZ_p& x, const vec_ZZ_p& a, long k,
                   const ZZ_pXArgument& H, const ZZ_pXModulus& F)

{
   long n = F.n;

   if (a.length() > n || k < 0 || k >= (1L << (NTL_BITS_PER_LONG-4))) 
      Error("ProjectPowers: bad args");

   long m = H.H.length()-1;
   long l = (k+m-1)/m - 1;

   ZZ_pXMultiplier M;
   build(M, H.H[m], F);

   vec_ZZ_p s(INIT_SIZE, n);
   s = a;
   StripZeroes(s);

   x.SetLength(k);

   for (long i = 0; i <= l; i++) {
      long m1 = min(m, k-i*m);
      ZZ_p* w = &x[i*m];
      for (long j = 0; j < m1; j++)
         InnerProduct(w[j], H.H[j].rep, s);
      if (i < l)
         UpdateMap(s, s, M, F);
   }
}



void ProjectPowers(vec_ZZ_p& x, const vec_ZZ_p& a, long k,
                   const ZZ_pX& h, const ZZ_pXModulus& F)

{
   if (a.length() > F.n || k < 0) Error("ProjectPowers: bad args");

   if (k == 0) {
      x.SetLength(0);
      return;
   }

   long m = SqrRoot(k);

   ZZ_pXArgument H;

   build(H, h, F, m);
   ProjectPowers(x, a, k, H, F);
}


void BerlekampMassey(ZZ_pX& h, const vec_ZZ_p& a, long m)
{
   ZZ_pX Lambda, Sigma, Temp;
   long L;
   ZZ_p Delta, Delta1, t1;
   long shamt;

   // cerr << "*** " << m << "\n";

   Lambda.SetMaxLength(m+1);
   Sigma.SetMaxLength(m+1);
   Temp.SetMaxLength(m+1);

   L = 0;
   set(Lambda);
   clear(Sigma);
   set(Delta);
   shamt = 0;

   long i, r, dl;

   for (r = 1; r <= 2*m; r++) {
      // cerr << r << "--";
      clear(Delta1);
      dl = deg(Lambda);
      for (i = 0; i <= dl; i++) {
         mul(t1, Lambda.rep[i], a[r-i-1]);
         add(Delta1, Delta1, t1);
      }

      if (IsZero(Delta1)) {
         shamt++;
         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else if (2*L < r) {
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         Sigma = Lambda;
         ShiftSub(Lambda, Temp, shamt+1);
         shamt = 0;
         L = r-L;
         Delta = Delta1;
         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else {
         shamt++;
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         ShiftSub(Lambda, Temp, shamt);
         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
   }

   // cerr << "finished: " << L << " " << deg(Lambda) << "\n"; 

   dl = deg(Lambda);
   h.rep.SetLength(L + 1);

   for (i = 0; i < L - dl; i++)
      clear(h.rep[i]);

   for (i = L - dl; i <= L; i++)
      h.rep[i] = Lambda.rep[L - i];
}




void GCDMinPolySeq(ZZ_pX& h, const vec_ZZ_p& x, long m)
{
   long i;
   ZZ_pX a, b;
   ZZ_pXMatrix M;
   ZZ_p t;

   a.rep.SetLength(2*m);
   for (i = 0; i < 2*m; i++) a.rep[i] = x[2*m-1-i];
   a.normalize();

   SetCoeff(b, 2*m);

   HalfGCD(M, b, a, m+1);

   /* make monic */

   inv(t, LeadCoeff(M(1,1)));
   mul(h, M(1,1), t);
}


void MinPolySeq(ZZ_pX& h, const vec_ZZ_p& a, long m)
{
   if (m < 0 || m >= (1L << (NTL_BITS_PER_LONG-4))) Error("MinPoly: bad args");
   if (a.length() < 2*m) Error("MinPoly: sequence too short");

   if (m > NTL_ZZ_pX_BERMASS_CROSSOVER)
      GCDMinPolySeq(h, a, m);
   else
      BerlekampMassey(h, a, m);
}


void DoMinPolyMod(ZZ_pX& h, const ZZ_pX& g, const ZZ_pXModulus& F, long m,
               const vec_ZZ_p& R) 
{
   vec_ZZ_p x;

   ProjectPowers(x, R, 2*m, g, F);
   MinPolySeq(h, x, m);
}