zz_px.c

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

      ZZ_pX P1;
      mul(P1, a, b);
      rem(x, P1, F);
      return;
   }

   d = da + db + 1;

   k = NextPowerOfTwo(d);
   k = max(k, F.k);

   FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, F.l);
   ZZ_pX P1(INIT_SIZE, n);

   ToFFTRep(R1, a, k);
   ToFFTRep(R2, b, k);
   mul(R1, R1, R2);
   NDFromFFTRep(P1, R1, n, d-1, R2); // save R1 for future use
   
   ToFFTRep(R2, P1, F.l);
   mul(R2, R2, F.HRep);
   FromFFTRep(P1, R2, n-2, 2*n-4);

   ToFFTRep(R2, P1, F.k);
   mul(R2, R2, F.FRep);
   reduce(R1, R1, F.k);
   sub(R1, R1, R2);
   FromFFTRep(x, R1, 0, n-1);
}

void SqrMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long  da, d, n, k;

   da = deg(a);
   n = F.n;

   if (n < 0) Error("SqrMod: uninitailized modulus");

   if (da >= n) 
      Error("bad args to SqrMod(ZZ_pX,ZZ_pX,ZZ_pXModulus)");

   if (!F.UseFFT || da <= NTL_ZZ_pX_FFT_CROSSOVER) {
      ZZ_pX P1;
      sqr(P1, a);
      rem(x, P1, F);
      return;
   }


   d = 2*da + 1;

   k = NextPowerOfTwo(d);
   k = max(k, F.k);

   FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, F.l);
   ZZ_pX P1(INIT_SIZE, n);

   ToFFTRep(R1, a, k);
   mul(R1, R1, R1);
   NDFromFFTRep(P1, R1, n, d-1, R2);  // save R1 for future use
   
   ToFFTRep(R2, P1, F.l);
   mul(R2, R2, F.HRep);
   FromFFTRep(P1, R2, n-2, 2*n-4);

   ToFFTRep(R2, P1, F.k);
   mul(R2, R2, F.FRep);
   reduce(R1, R1, F.k);
   sub(R1, R1, R2);
   FromFFTRep(x, R1, 0, n-1);
}

void PlainInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m)

   /* x = (1/a) % X^m, input not output, constant term a is nonzero */

{
   long i, k, n, lb;
   static ZZ v, t;
   ZZ_p s;
   const ZZ_p* ap;
   ZZ_p* xp;
   

   n = deg(a);

   if (n < 0) Error("division by zero");

   inv(s, ConstTerm(a));

   if (n == 0) {
      conv(x, s);
      return;
   }

   ap = a.rep.elts();
   x.rep.SetLength(m);
   xp = x.rep.elts();

   xp[0] = s;

   long is_one = IsOne(s);

   for (k = 1; k < m; k++) {
      clear(v);
      lb = max(k-n, 0);
      for (i = lb; i <= k-1; i++) {
         mul(t, rep(xp[i]), rep(ap[k-i]));
         add(v, v, t);
      }
      conv(xp[k], v);
      negate(xp[k], xp[k]);
      if (!is_one) mul(xp[k], xp[k], s);
   }

   x.normalize();
}


void trunc(ZZ_pX& x, const ZZ_pX& a, long m)

// x = a % X^m, output may alias input 

{
   if (m < 0) Error("trunc: bad args");

   if (&x == &a) {
      if (x.rep.length() > m) {
         x.rep.SetLength(m);
         x.normalize();
      }
   }
   else {
      long n;
      long i;
      ZZ_p* xp;
      const ZZ_p* ap;

      n = min(a.rep.length(), m);
      x.rep.SetLength(n);

      xp = x.rep.elts();
      ap = a.rep.elts();

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

      x.normalize();
   }
}

void CyclicReduce(ZZ_pX& x, const ZZ_pX& a, long m)

// computes x = a mod X^m-1

{
   long n = deg(a);
   long i, j;
   ZZ_p accum;

   if (n < m) {
      x = a;
      return;
   }

   if (&x != &a)
      x.rep.SetLength(m);

   for (i = 0; i < m; i++) {
      accum = a.rep[i];
      for (j = i + m; j <= n; j += m)
         add(accum, accum, a.rep[j]);
      x.rep[i] = accum;
   }

   if (&x == &a)
      x.rep.SetLength(m);

   x.normalize();
}



void InvTrunc(ZZ_pX& x, const ZZ_pX& a, long m)
{
   if (m < 0) Error("InvTrunc: bad args");

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

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

   if (&x == &a) {
      ZZ_pX la;
      la = a;
      if (m > NTL_ZZ_pX_NEWTON_CROSSOVER && deg(a) > 0)
         NewtonInvTrunc(x, la, m);
      else
         PlainInvTrunc(x, la, m);
   }
   else {
      if (m > NTL_ZZ_pX_NEWTON_CROSSOVER && deg(a) > 0)
         NewtonInvTrunc(x, a, m);
      else
         PlainInvTrunc(x, a, m);
   }
}
   


void build(ZZ_pXModulus& x, const ZZ_pX& f)
{
   x.f = f;
   x.n = deg(f);

   x.tracevec.SetLength(0);

   if (x.n <= 0)
      Error("build: deg(f) must be at least 1");

   if (x.n <= NTL_ZZ_pX_FFT_CROSSOVER + 1) {
      x.UseFFT = 0;
      return;
   }

   x.UseFFT = 1;

   x.k = NextPowerOfTwo(x.n);
   x.l = NextPowerOfTwo(2*x.n - 3);
   ToFFTRep(x.FRep, f, x.k);

   ZZ_pX P1(INIT_SIZE, x.n+1), P2(INIT_SIZE, x.n);

   CopyReverse(P1, f, 0, x.n);
   InvTrunc(P2, P1, x.n-1);

   CopyReverse(P1, P2, 0, x.n-2);
   ToFFTRep(x.HRep, P1, x.l);
}

ZZ_pXModulus::ZZ_pXModulus(const ZZ_pX& ff)
{
   build(*this, ff);
}

ZZ_pXMultiplier::ZZ_pXMultiplier(const ZZ_pX& b, const ZZ_pXModulus& F)
{
   build(*this, b, F); 
}

void build(ZZ_pXMultiplier& x, const ZZ_pX& b, 
                         const ZZ_pXModulus& F)
{
   long db;
   long n = F.n;

   if (n < 0) Error("build ZZ_pXMultiplier: uninitialized modulus");

   x.b = b;
   db = deg(b);

   if (db >= n) Error("build ZZ_pXMultiplier: deg(b) >= deg(f)");

   if (!F.UseFFT || db <= NTL_ZZ_pX_FFT_CROSSOVER) {
      x.UseFFT = 0;
      return;
   }

   x.UseFFT = 1;

   FFTRep R1(INIT_SIZE, F.l);
   ZZ_pX P1(INIT_SIZE, n);
   

   ToFFTRep(R1, b, F.l);
   reduce(x.B2, R1, F.k);
   mul(R1, R1, F.HRep);
   FromFFTRep(P1, R1, n-1, 2*n-3); 
   ToFFTRep(x.B1, P1, F.l);
}


void MulMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXMultiplier& B,
                                      const ZZ_pXModulus& F)
{

   long n = F.n;
   long da;

   da = deg(a);

   if (da >= n)
      Error(" bad args to MulMod(ZZ_pX,ZZ_pX,ZZ_pXMultiplier,ZZ_pXModulus)");

   if (da < 0) {
      clear(x);
      return;
   }

   if (!B.UseFFT || !F.UseFFT || da <= NTL_ZZ_pX_FFT_CROSSOVER) {
      ZZ_pX P1;
      mul(P1, a, B.b);
      rem(x, P1, F);
      return;
   }

   ZZ_pX P1(INIT_SIZE, n), P2(INIT_SIZE, n);
   FFTRep R1(INIT_SIZE, F.l), R2(INIT_SIZE, F.l);

   ToFFTRep(R1, a, F.l);
   mul(R2, R1, B.B1);
   FromFFTRep(P1, R2, n-1, 2*n-3);

   reduce(R1, R1, F.k);
   mul(R1, R1, B.B2);
   ToFFTRep(R2, P1, F.k);
   mul(R2, R2, F.FRep);
   sub(R1, R1, R2);

   FromFFTRep(x, R1, 0, n-1);
}
   

void PowerXMod(ZZ_pX& hh, const ZZ& e, const ZZ_pXModulus& F)
{
   if (F.n < 0) Error("PowerXMod: uninitialized modulus");

   if (IsZero(e)) {
      set(hh);
      return;
   }

   long n = NumBits(e);
   long i;

   ZZ_pX h;

   h.SetMaxLength(F.n);
   set(h);

   for (i = n - 1; i >= 0; i--) {
      SqrMod(h, h, F);
      if (bit(e, i))
         MulByXMod(h, h, F);
   }

   if (e < 0) InvMod(h, h, F);

   hh = h;
}


void PowerXPlusAMod(ZZ_pX& hh, const ZZ_p& a, const ZZ& e, const ZZ_pXModulus& F)
{
   if (F.n < 0) Error("PowerXPlusAMod: uninitialized modulus");

   if (IsZero(e)) {
      set(hh);
      return;
   }

   ZZ_pX t1(INIT_SIZE, F.n), t2(INIT_SIZE, F.n);
   long n = NumBits(e);
   long i;

   ZZ_pX h;

   h.SetMaxLength(F.n);
   set(h);

   for (i = n - 1; i >= 0; i--) {
      SqrMod(h, h, F);
      if (bit(e, i)) {
         MulByXMod(t1, h, F);
         mul(t2, h, a);
         add(h, t1, t2);
      }
   }

   if (e < 0) InvMod(h, h, F);

   hh = h;
}


void PowerMod(ZZ_pX& h, const ZZ_pX& g, const ZZ& e, const ZZ_pXModulus& F)
{
   if (deg(g) >= F.n)
      Error("PowerMod: bad args");

   if (IsZero(e)) {
      set(h);
      return;
   }

   ZZ_pXMultiplier G;

   ZZ_pX res;

   long n = NumBits(e);
   long i;

   build(G, g, F);

   res.SetMaxLength(F.n);
   set(res);

   for (i = n - 1; i >= 0; i--) {
      SqrMod(res, res, F);
      if (bit(e, i))
         MulMod(res, res, G, F);
   }

   if (e < 0) InvMod(res, res, F);

   h = res;
}

#if 0
void NewtonInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m)
{
   x.SetMaxLength(m);

   long i;
   long t;


   t = NextPowerOfTwo(2*m-1);

   FFTRep R1(INIT_SIZE, t), R2(INIT_SIZE, t);
   ZZ_pX P1(INIT_SIZE, m);


   long log2_newton = NextPowerOfTwo(NTL_ZZ_pX_NEWTON_CROSSOVER)-1;
   PlainInv(x, a, 1L << log2_newton);
   long k = 1L << log2_newton;
   long a_len = min(m, a.rep.length());

   while (k < m) {
      long l = min(2*k, m);

      t = NextPowerOfTwo(2*k);
      ToFFTRep(R1, x, t);
      mul(R1, R1, R1);
      FromFFTRep(P1, R1, 0, l-1);

      t = NextPowerOfTwo(deg(P1) + min(l, a_len));
      ToFFTRep(R1, P1, t);
      ToFFTRep(R2, a, t, 0, min(l, a_len)-1);
      mul(R1, R1, R2);
      FromFFTRep(P1, R1, k, l-1);
      
      x.rep.SetLength(l);
      long y_len = P1.rep.length();
      for (i = k; i < l; i++) {
         if (i-k >= y_len)
            clear(x.rep[i]);
         else
            negate(x.rep[i], P1.rep[i-k]);
      }
      x.normalize();

      k = l;
   }
}


#else

void NewtonInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m)
{
   x.SetMaxLength(m);
   long i, t, k;

   long log2_newton = NextPowerOfTwo(NTL_ZZ_pX_NEWTON_CROSSOVER)-1;
   PlainInvTrunc(x, a, 1L << log2_newton);

   t = NextPowerOfTwo(m);

   FFTRep R1(INIT_SIZE, t), R2(INIT_SIZE, t);
   ZZ_pX P1(INIT_SIZE, m/2);

   long a_len = min(m, a.rep.length());

   ZZ_pXModRep a_rep;
   ToZZ_pXModRep(a_rep, a, 0, a_len-1);

   k = 1L << log2_newton; 
   t = log2_newton;

   while (k < m) {
      long l = min(2*k, m);

      ToFFTRep(R1, x, t+1);
      ToFFTRep(R2, a_rep, t+1, 0, l-1); 
      mul(R2, R2, R1);
      FromFFTRep(P1, R2, k, l-1);
      
      ToFFTRep(R2, P1, t+1);
      mul(R2, R2, R1);
      FromFFTRep(P1, R2, 0, l-k-1);

      x.rep.SetLength(l);
      long y_len = P1.rep.length();
      for (i = k; i < l; i++) {
         if (i-k >= y_len)
            clear(x.rep[i]);
         else
            negate(x.rep[i], P1.rep[i-k]);
      }
      x.normalize();

      t++;
      k = l;
   }
}



#endif


void FFTDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b)
{
   long n = deg(b);
   long m = deg(a);
   long k, l;

   if (m < n) {
      clear(q);
      r = a;
      return;
   }

   if (m >= 3*n) {
      ZZ_pXModulus B;
      build(B, b);
      DivRem(q, r, a, B);
      return;
   }

   ZZ_pX P1, P2, P3;

   CopyReverse(P3, b, 0, n);
   InvTrunc(P2, P3, m-n+1);
   CopyReverse(P1, P2, 0, m-n);

   k = NextPowerOfTwo(2*(m-n)+1);
   long k1 = NextPowerOfTwo(n);
   long mx = max(k1, k);

   FFTRep R1(INIT_SIZE, mx), R2(INIT_SIZE, mx);

   ToFFTRep(R1, P1, k);
   ToFFTRep(R2, a, k, n, m);
   mul(R1, R1, R2);
   FromFFTRep(P3, R1, m-n, 2*(m-n));
   
   l = 1L << k1;

   
   ToFFTRep(R1, b, k1);
   ToFFTRep(R2, P3, k1);
   mul(R1, R1, R2);
   FromFFTRep(P1, R1, 0, n-1);
   CyclicReduce(P2, a, l);
   trunc(r, P2, n);
   sub(r, r, P1);
   q = P3;
}




void FFTDiv(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b)
{

   long n = deg(b);
   long m = deg(a);
   long k;

   if (m < n) {
      clear(q);
      return;
   }

   if (m >= 3*n) {
      ZZ_pXModulus B;
      build(B, b);
      div(q, a, B);
      return;
   }

   ZZ_pX P1, P2, P3;

   CopyReverse(P3, b, 0, n);
   InvTrunc(P2, P3, m-n+1);
   CopyReverse(P1, P2, 0, m-n);

   k = NextPowerOfTwo(2*(m-n)+1);

   FFTRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);

   ToFFTRep(R1, P1, k);
   ToFFTRep(R2, a, k, n, m);
   mul(R1, R1, R2);
   FromFFTRep(q, R1, m-n, 2*(m-n));
}



void FFTRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b)
{
   long n = deg(b);
   long m = deg(a);
   long k, l;

   if (m < n) {
      r = a;
      return;
   }

   if (m >= 3*n) {
      ZZ_pXModulus B;
      build(B, b);
      rem(r, a, B);
      return;
   }

   ZZ_pX P1, P2, P3;

   CopyReverse(P3, b, 0, n);
   InvTrunc(P2, P3, m-n+1);
   CopyReverse(P1, P2, 0, m-n);

   k = NextPowerOfTwo(2*(m-n)+1);
   long k1 = NextPowerOfTwo(n);
   long mx = max(k, k1);

   FFTRep R1(INIT_SIZE, mx), R2(INIT_SIZE, mx);

   ToFFTRep(R1, P1, k);
   ToFFTRep(R2, a, k, n, m);
   mul(R1, R1, R2);
   FromFFTRep(P3, R1, m-n, 2*(m-n));
   
   l = 1L << k1;
   
   ToFFTRep(R1, b, k1);
   ToFFTRep(R2, P3, k1);
   mul(R1, R1, R2);
   FromFFTRep(P3, R1, 0, n-1);
   CyclicReduce(P2, a, l);
   trunc(r, P2, n);
   sub(r, r, P3);
}


void DivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b)
{
   if (deg(b) > NTL_ZZ_pX_DIV_CROSSOVER && deg(a) - deg(b) > NTL_ZZ_pX_DIV_CROSSOVER)
      FFTDivRem(q, r, a, b);
   else
      PlainDivRem(q, r, a, b);
}

void div(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b)
{
   if (deg(b) > NTL_ZZ_pX_DIV_CROSSOVER && deg(a) - deg(b) > NTL_ZZ_pX_DIV_CROSSOVER)
      FFTDiv(q, a, b);
   else
      PlainDiv(q, a, b);
}

void div(ZZ_pX& q, const ZZ_pX& a, const ZZ_p& b)
{
   ZZ_pTemp TT; ZZ_p& T = TT.val();

   inv(T, b);
   mul(q, a, T);
}

void div(ZZ_pX& q, const ZZ_pX& a, long b)
{
   ZZ_pTemp TT; ZZ_p& T = TT.val();

   T = b;
   inv(T, T);
   mul(q, a, T);
}



void rem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b)
{
   if (deg(b) > NTL_ZZ_pX_DIV_CROSSOVER && deg(a) - deg(b) > NTL_ZZ_pX_DIV_CROSSOVER)
      FFTRem(r, a, b);
   else
      PlainRem(r, a, b);
}


long operator==(const ZZ_pX& a, long b)
{
   if (b == 0)
      return IsZero(a);

   if (b == 1)
      return IsOne(a);

   long da = deg(a);

   if (da > 0)
      return 0;

   ZZ_pTemp TT; ZZ_p& bb = TT.val();
   bb = b;

   if (da < 0)
      return IsZero(bb);

   return a.rep[0] == bb;
}

long operator==(const ZZ_pX& a, const ZZ_p& b)
{
   if (IsZero(b))
      return IsZero(a);

   long da = deg(a);

   if (da != 0)
      return 0;

   return a.rep[0] == b;
}

void power(ZZ_pX& x, const ZZ_pX& a, long e)
{
   if (e < 0) {
      Error("power: negative exponent");
   }

   if (e == 0) {
      x = 1;
      return;
   }

   if (a == 0 || a == 1) {
      x = a;
      return;
   }

   long da = deg(a);

   if (da == 0) {
      x = power(ConstTerm(a), e);
      return;
   }

   if (da > (NTL_MAX_LONG-1)/e)
      Error("overflow in power");

   ZZ_pX res;
   res.SetMaxLength(da*e + 1);
   res = 1;
   
   long k = NumBits(e);
   long i;

   for (i = k - 1; i >= 0; i--) {
      sqr(res, res);
      if (bit(e, i))
         mul(res, res, a);
   }

   x = res;
}

void reverse(ZZ_pX& x, const ZZ_pX& a, long hi)
{
   if (hi < -1) Error("reverse: bad args");
   if (hi >= (1L << (NTL_BITS_PER_LONG-4)))
      Error("overflow in reverse");

   if (&x == &a) {
      ZZ_pX tmp;
      CopyReverse(tmp, a, 0, hi);
      x = tmp;
   }
   else
      CopyReverse(x, a, 0, hi);
}

NTL_END_IMPL