zz_px.c

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

}

void FromFFTRep(ZZ_pX& x, FFTRep& y, long lo, long hi)

   // converts from FFT-representation to coefficient representation
   // only the coefficients lo..hi are computed
   

{
   ZZ_pInfo->check();

   long k, n, i, j, l;

   vec_long& t = ModularRepBuf;
   vec_long& s = FFTBuf;;

   t.SetLength(ZZ_pInfo->NumPrimes);

   k = y.k;
   n = (1L << k);

   s.SetLength(n);
   long *sp = s.elts();

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *yp = &y.tbl[i][0];
      long q = FFTPrime[i];
      double qinv = ((double) 1)/((double) q);
      long w = TwoInvTable[i][k];
      long *Root = &RootInvTable[i][0];

      FFT(sp, yp, k, q, Root);

      for (j = 0; j < n; j++) yp[j] = MulMod(sp[j], w, q, qinv);
   }

   hi = min(hi, n-1);
   l = hi-lo+1;
   l = max(l, 0);
   x.rep.SetLength(l);

   for (j = 0; j < l; j++) {
      for (i = 0; i < ZZ_pInfo->NumPrimes; i++) 
         t[i] = y.tbl[i][j+lo]; 

      FromModularRep(x.rep[j], t);
   }

   x.normalize();
}

void RevFromFFTRep(vec_ZZ_p& x, FFTRep& y, long lo, long hi)

   // converts from FFT-representation to coefficient representation
   // using "inverted" evaluation points.
   // only the coefficients lo..hi are computed
   

{
   ZZ_pInfo->check();

   long k, n, i, j, l;

   vec_long& t = ModularRepBuf;
   vec_long& s = FFTBuf;

   k = y.k;
   n = (1L << k);

   t.SetLength(ZZ_pInfo->NumPrimes);
   s.SetLength(n);
   long *sp = s.elts();

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *yp = &y.tbl[i][0];
      long q = FFTPrime[i];
      long *Root = &RootTable[i][0];

      FFT(sp, yp, k, q, Root);
      for (j = 0; j < n; j++)
         yp[j] = sp[j];
   }

   hi = min(hi, n-1);
   l = hi-lo+1;
   l = max(l, 0);
   x.SetLength(l);

   for (j = 0; j < l; j++) {
      for (i = 0; i < ZZ_pInfo->NumPrimes; i++) 
         t[i] = y.tbl[i][j+lo]; 

      FromModularRep(x[j], t);
   }
}

void NDFromFFTRep(ZZ_pX& x, const FFTRep& y, long lo, long hi, FFTRep& z)
{
   ZZ_pInfo->check();

   long k, n, i, j, l;

   vec_long& t = ModularRepBuf;

   t.SetLength(ZZ_pInfo->NumPrimes);
   k = y.k;
   n = (1L << k);

   z.SetSize(k);

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *zp = &z.tbl[i][0];
      long q = FFTPrime[i];
      double qinv = ((double) 1)/((double) q);
      long w = TwoInvTable[i][k];
      long *Root = &RootInvTable[i][0];

      FFT(zp, &y.tbl[i][0], k, q, Root);

      for (j = 0; j < n; j++) zp[j] = MulMod(zp[j], w, q, qinv);
   }

   hi = min(hi, n-1);
   l = hi-lo+1;
   l = max(l, 0);
   x.rep.SetLength(l);

   for (j = 0; j < l; j++) {
      for (i = 0; i < ZZ_pInfo->NumPrimes; i++) 
         t[i] = z.tbl[i][j+lo]; 

      FromModularRep(x.rep[j], t);
   }

   x.normalize();
}

void NDFromFFTRep(ZZ_pX& x, FFTRep& y, long lo, long hi)
{
   FFTRep z;
   NDFromFFTRep(x, y, lo, hi, z);
}

void FromFFTRep(ZZ_p* x, FFTRep& y, long lo, long hi)

   // converts from FFT-representation to coefficient representation
   // only the coefficients lo..hi are computed
   

{
   ZZ_pInfo->check();

   long k, n, i, j;

   vec_long& t = ModularRepBuf;
   vec_long& s = FFTBuf;

   k = y.k;
   n = (1L << k);

   t.SetLength(ZZ_pInfo->NumPrimes);
   s.SetLength(n);
   long *sp = s.elts();

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *yp = &y.tbl[i][0];
      long q = FFTPrime[i];
      double qinv = ((double) 1)/((double) q);
      long w = TwoInvTable[i][k];
      long *Root = &RootInvTable[i][0];

      FFT(sp, yp, k, q, Root);

      for (j = 0; j < n; j++) yp[j] = MulMod(sp[j], w, q, qinv);
   }

   for (j = lo; j <= hi; j++) {
      if (j >= n)
         clear(x[j-lo]);
      else {
         for (i = 0; i < ZZ_pInfo->NumPrimes; i++) 
            t[i] = y.tbl[i][j]; 

         FromModularRep(x[j-lo], t);
      }
   }
}


void mul(FFTRep& z, const FFTRep& x, const FFTRep& y)
{
   ZZ_pInfo->check();

   long k, n, i, j;

   if (x.k != y.k) Error("FFT rep mismatch");

   k = x.k;
   n = 1L << k;

   z.SetSize(k);

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *zp = &z.tbl[i][0];
      const long *xp = &x.tbl[i][0];
      const long *yp = &y.tbl[i][0];
      long q = FFTPrime[i];
      double qinv = ((double) 1)/((double) q);

      for (j = 0; j < n; j++)
         zp[j] = MulMod(xp[j], yp[j], q, qinv);
   }

}

void sub(FFTRep& z, const FFTRep& x, const FFTRep& y)
{
   ZZ_pInfo->check();

   long k, n, i, j;

   if (x.k != y.k) Error("FFT rep mismatch");

   k = x.k;
   n = 1L << k;

   z.SetSize(k);

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *zp = &z.tbl[i][0];
      const long *xp = &x.tbl[i][0];
      const long *yp = &y.tbl[i][0];
      long q = FFTPrime[i];

      for (j = 0; j < n; j++)
         zp[j] = SubMod(xp[j], yp[j], q);
   }
}

void add(FFTRep& z, const FFTRep& x, const FFTRep& y)
{
   ZZ_pInfo->check();

   long k, n, i, j;

   if (x.k != y.k) Error("FFT rep mismatch");

   k = x.k;
   n = 1L << k;

   z.SetSize(k);

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long *zp = &z.tbl[i][0];
      const long *xp = &x.tbl[i][0];
      const long *yp = &y.tbl[i][0];
      long q = FFTPrime[i];

      for (j = 0; j < n; j++)
         zp[j] = AddMod(xp[j], yp[j], q);
   }
}


void reduce(FFTRep& x, const FFTRep& a, long k)
  // reduces a 2^l point FFT-rep to a 2^k point FFT-rep
  // input may alias output
{
   ZZ_pInfo->check();

   long i, j, l, n;
   long* xp;
   const long* ap;

   l = a.k;
   n = 1L << k;

   if (l < k) Error("reduce: bad operands");

   x.SetSize(k);

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      ap = &a.tbl[i][0];   
      xp = &x.tbl[i][0];
      for (j = 0; j < n; j++) 
         xp[j] = ap[j << (l-k)];
   }
}

void AddExpand(FFTRep& x, const FFTRep& a)
//  x = x + (an "expanded" version of a)
{
   ZZ_pInfo->check();

   long i, j, l, k, n;

   l = x.k;
   k = a.k;
   n = 1L << k;

   if (l < k) Error("AddExpand: bad args");

   for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
      long q = FFTPrime[i];
      const long *ap = &a.tbl[i][0];
      long *xp = &x.tbl[i][0];
      for (j = 0; j < n; j++) {
         long j1 = j << (l-k);
         xp[j1] = AddMod(xp[j1], ap[j], q);
      }
   }
}



void ToZZ_pXModRep(ZZ_pXModRep& y, const ZZ_pX& x, long lo, long hi)
{
   ZZ_pInfo->check();

   long n, i, j;
   vec_long& t = ModularRepBuf;

   t.SetLength(ZZ_pInfo->NumPrimes);

   if (lo < 0)
      Error("bad arg to ToZZ_pXModRep");
   hi = min(hi, deg(x));
   n = max(hi-lo+1, 0);

   y.SetSize(n);

   const ZZ_p *xx = x.rep.elts();

   for (j = 0; j < n; j++) {
      ToModularRep(t, xx[j+lo]);
      for (i = 0; i < ZZ_pInfo->NumPrimes; i++)
         y.tbl[i][j] = t[i];
   }
}


void ToFFTRep(FFTRep& x, const ZZ_pXModRep& a, long k, long lo, long hi)
{
   ZZ_pInfo->check();

   vec_long s;
   long n, m, i, j;

   if (k < 0 || lo < 0)
      Error("bad args to ToFFTRep");

   if (hi > a.n-1) hi = a.n-1;

   n = 1L << k;
   m = max(hi-lo+1, 0);

   if (m > n)
      Error("bad args to ToFFTRep");

   s.SetLength(n);
   long *sp = s.elts();

   x.SetSize(k);

   long NumPrimes = ZZ_pInfo->NumPrimes;

   for (i = 0; i < NumPrimes; i++) {
      long *Root = &RootTable[i][0];
      long *xp = &x.tbl[i][0];
      long *ap = (m == 0 ? 0 : &a.tbl[i][0]);
      for (j = 0; j < m; j++)
         sp[j] = ap[lo+j];
      for (j = m; j < n; j++)
         sp[j] = 0;
      
      FFT(xp, sp, k, FFTPrime[i], Root);
   }
}






void FFTMul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b)
{
   long k, d;

   if (IsZero(a) || IsZero(b)) {
      clear(x);
      return;
   }

   d = deg(a) + deg(b);
   k = NextPowerOfTwo(d+1);

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

   ToFFTRep(R1, a, k);
   ToFFTRep(R2, b, k);
   mul(R1, R1, R2);
   FromFFTRep(x, R1, 0, d);
}

void FFTSqr(ZZ_pX& x, const ZZ_pX& a)
{
   long k, d;

   if (IsZero(a)) {
      clear(x);
      return;
   }

   d = 2*deg(a);
   k = NextPowerOfTwo(d+1);

   FFTRep R1(INIT_SIZE, k);

   ToFFTRep(R1, a, k);
   mul(R1, R1, R1);
   FromFFTRep(x, R1, 0, d);
}


void CopyReverse(ZZ_pX& x, const ZZ_pX& a, long lo, long hi)

   // x[0..hi-lo] = reverse(a[lo..hi]), with zero fill
   // input may not alias output

{
   long i, j, n, m;

   n = hi-lo+1;
   m = a.rep.length();

   x.rep.SetLength(n);

   const ZZ_p* ap = a.rep.elts();
   ZZ_p* xp = x.rep.elts();

   for (i = 0; i < n; i++) {
      j = hi-i;
      if (j < 0 || j >= m)
         clear(xp[i]);
      else
         xp[i] = ap[j];
   }

   x.normalize();
} 

void copy(ZZ_pX& x, const ZZ_pX& a, long lo, long hi)

   // x[0..hi-lo] = a[lo..hi], with zero fill
   // input may not alias output

{
   long i, j, n, m;

   n = hi-lo+1;
   m = a.rep.length();

   x.rep.SetLength(n);

   const ZZ_p* ap = a.rep.elts();
   ZZ_p* xp = x.rep.elts();

   for (i = 0; i < n; i++) {
      j = lo + i;
      if (j < 0 || j >= m)
         clear(xp[i]);
      else
         xp[i] = ap[j];
   }

   x.normalize();
} 


void rem21(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long i, da, ds, n, kk;

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

   if (da > 2*n-2)
      Error("bad args to rem(ZZ_pX,ZZ_pX,ZZ_pXModulus)");


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

   if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) {
      PlainRem(x, a, F.f);
      return;
   }

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

   ToFFTRep(R1, a, F.l, n, 2*(n-1));
   mul(R1, R1, F.HRep);
   FromFFTRep(P1, R1, n-2, 2*n-4);

   ToFFTRep(R1, P1, F.k);
   mul(R1, R1, F.FRep);
   FromFFTRep(P1, R1, 0, n-1);

   ds = deg(P1);

   kk = 1L << F.k;

   x.rep.SetLength(n);
   const ZZ_p* aa = a.rep.elts();
   const ZZ_p* ss = P1.rep.elts();
   ZZ_p* xx = x.rep.elts();

   for (i = 0; i < n; i++) {
      if (i <= ds)
         sub(xx[i], aa[i], ss[i]);
      else
         xx[i] = aa[i];

      if (i + kk <= da)
         add(xx[i], xx[i], aa[i+kk]);
   }

   x.normalize();
}

void DivRem21(ZZ_pX& q, ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long i, da, ds, n, kk;

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

   if (da > 2*n-2)
      Error("bad args to rem(ZZ_pX,ZZ_pX,ZZ_pXModulus)");


   if (da < n) {
      x = a;
      clear(q);
      return;
   }

   if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) {
      PlainDivRem(q, x, a, F.f);
      return;
   }

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

   ToFFTRep(R1, a, F.l, n, 2*(n-1));
   mul(R1, R1, F.HRep);
   FromFFTRep(P1, R1, n-2, 2*n-4);
   qq = P1;

   ToFFTRep(R1, P1, F.k);
   mul(R1, R1, F.FRep);
   FromFFTRep(P1, R1, 0, n-1);

   ds = deg(P1);

   kk = 1L << F.k;

   x.rep.SetLength(n);
   const ZZ_p* aa = a.rep.elts();
   const ZZ_p* ss = P1.rep.elts();
   ZZ_p* xx = x.rep.elts();

   for (i = 0; i < n; i++) {
      if (i <= ds)
         sub(xx[i], aa[i], ss[i]);
      else
         xx[i] = aa[i];

      if (i + kk <= da)
         add(xx[i], xx[i], aa[i+kk]);
   }

   x.normalize();
   q = qq;
}

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

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

   if (da > 2*n-2)
      Error("bad args to rem(ZZ_pX,ZZ_pX,ZZ_pXModulus)");


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

   if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) {
      PlainDiv(x, a, F.f);
      return;
   }

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

   ToFFTRep(R1, a, F.l, n, 2*(n-1));
   mul(R1, R1, F.HRep);
   FromFFTRep(x, R1, n-2, 2*n-4);
}


void rem(ZZ_pX& x, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long da = deg(a);
   long n = F.n;

   if (n < 0) Error("rem: unitialized modulus");

   if (da <= 2*n-2) {
      rem21(x, a, F);
      return;
   }
   else if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) {
      PlainRem(x, a, F.f);
      return;
   }

   ZZ_pX buf(INIT_SIZE, 2*n-1);

   long a_len = da+1;

   while (a_len > 0) {
      long old_buf_len = buf.rep.length();
      long amt = min(2*n-1-old_buf_len, a_len);

      buf.rep.SetLength(old_buf_len+amt);

      long i;

      for (i = old_buf_len+amt-1; i >= amt; i--)
         buf.rep[i] = buf.rep[i-amt];

      for (i = amt-1; i >= 0; i--)
         buf.rep[i] = a.rep[a_len-amt+i];

      buf.normalize();

      rem21(buf, buf, F);

      a_len -= amt;
   }

   x = buf;
}

void DivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long da = deg(a);
   long n = F.n;

   if (n < 0) Error("uninitialized modulus");

   if (da <= 2*n-2) {
      DivRem21(q, r, a, F);
      return;
   }
   else if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) {
      PlainDivRem(q, r, a, F.f);
      return;
   }

   ZZ_pX buf(INIT_SIZE, 2*n-1);
   ZZ_pX qbuf(INIT_SIZE, n-1);

   ZZ_pX qq;
   qq.rep.SetLength(da-n+1);

   long a_len = da+1;
   long q_hi = da-n+1;

   while (a_len > 0) {
      long old_buf_len = buf.rep.length();
      long amt = min(2*n-1-old_buf_len, a_len);

      buf.rep.SetLength(old_buf_len+amt);

      long i;

      for (i = old_buf_len+amt-1; i >= amt; i--)
         buf.rep[i] = buf.rep[i-amt];

      for (i = amt-1; i >= 0; i--)
         buf.rep[i] = a.rep[a_len-amt+i];

      buf.normalize();

      DivRem21(qbuf, buf, buf, F);
      long dl = qbuf.rep.length();
      a_len = a_len - amt;
      for(i = 0; i < dl; i++)
         qq.rep[a_len+i] = qbuf.rep[i];
      for(i = dl+a_len; i < q_hi; i++)
         clear(qq.rep[i]);
      q_hi = a_len;
   }

   r = buf;

   qq.normalize();
   q = qq;
}

void div(ZZ_pX& q, const ZZ_pX& a, const ZZ_pXModulus& F)
{
   long da = deg(a);
   long n = F.n;

   if (n < 0) Error("uninitialized modulus");

   if (da <= 2*n-2) {
      div21(q, a, F);
      return;
   }
   else if (!F.UseFFT || da - n <= NTL_ZZ_pX_FFT_CROSSOVER) {
      PlainDiv(q, a, F.f);
      return;
   }

   ZZ_pX buf(INIT_SIZE, 2*n-1);
   ZZ_pX qbuf(INIT_SIZE, n-1);

   ZZ_pX qq;
   qq.rep.SetLength(da-n+1);

   long a_len = da+1;
   long q_hi = da-n+1;

   while (a_len > 0) {
      long old_buf_len = buf.rep.length();
      long amt = min(2*n-1-old_buf_len, a_len);

      buf.rep.SetLength(old_buf_len+amt);

      long i;

      for (i = old_buf_len+amt-1; i >= amt; i--)
         buf.rep[i] = buf.rep[i-amt];

      for (i = amt-1; i >= 0; i--)
         buf.rep[i] = a.rep[a_len-amt+i];

      buf.normalize();

      a_len = a_len - amt;
      if (a_len > 0)
         DivRem21(qbuf, buf, buf, F);
      else
         div21(qbuf, buf, F);

      long dl = qbuf.rep.length();
      for(i = 0; i < dl; i++)
         qq.rep[a_len+i] = qbuf.rep[i];
      for(i = dl+a_len; i < q_hi; i++)
         clear(qq.rep[i]);
      q_hi = a_len;
   }

   qq.normalize();
   q = qq;
}



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

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

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

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

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

   if (!F.UseFFT || da <= NTL_ZZ_pX_FFT_CROSSOVER || db <= NTL_ZZ_pX_FFT_CROSSOVER) {