lzz_px.c

数值算法库for Unix

            y.tbl[0][offset] = 0;
         }
         else {
            accum = xx[j+lo];
            for (j1 = j + n; j1 < m; j1 += n)
               add(accum, accum, xx[j1+lo]);
               y.tbl[0][offset] = rep(accum);
         }
         offset = (offset + 1) & (n-1);
      }
   }
   else {
      for (j = 0; j < n; j++) {
         if (j >= m) {
            for (i = 0; i < NumPrimes; i++)
               y.tbl[i][offset] = 0;
         }
         else {
            accum = xx[j+lo];
            for (j1 = j + n; j1 < m; j1 += n)
               add(accum, accum, xx[j1+lo]);
            for (i = 0; i < NumPrimes; i++) {
               long q = FFTPrime[i];
               long t = rep(accum);
               if (t >= q) t -= q;
               y.tbl[i][offset] = t;
            }
         }
         offset = (offset + 1) & (n-1);
      }
   }


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

   if (index >= 0) {
      long *Root = &RootInvTable[index][0];
      long *yp = &y.tbl[0][0];
      long w = TwoInvTable[index][k];
      long q = FFTPrime[index];
      double qinv = ((double) 1)/((double) q);
      FFT(sp, yp, y.k, q, Root);
      for (j = 0; j < n; j++)
         yp[j] = MulMod(sp[j], w, q, qinv);
   }
   else {
      for (i = 0; i < zz_pInfo->NumPrimes; i++) {
         long *Root = &RootInvTable[i][0];
         long *yp = &y.tbl[i][0];
         long w = TwoInvTable[i][k];
         long q = FFTPrime[i];
         double qinv = ((double) 1)/((double) q);
         FFT(sp, yp, y.k, q, Root);
         for (j = 0; j < n; j++)
            yp[j] = MulMod(sp[j], w, q, qinv);
      }
   }
}

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
   

{
   long k, n, i, j, l;
   long NumPrimes = zz_pInfo->NumPrimes;

   long t[4];
   vec_long& s = FFTBuf;;

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

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

   long index = zz_pInfo->index;

   if (index >= 0) {
      long *yp = &y.tbl[0][0];
      long q = FFTPrime[index];
      double qinv = FFTPrimeInv[index];
      long w = TwoInvTable[index][k];
      long *Root = &RootInvTable[index][0];

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

      for (j = 0; j < n; j++) yp[j] = MulMod(sp[j], w, q, qinv);
   }
   else {
      for (i = 0; i < NumPrimes; i++) {
         long *yp = &y.tbl[i][0];
         long q = FFTPrime[i];
         double qinv = FFTPrimeInv[i];
         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);

   if (index >= 0) {
      zz_p *xp = x.rep.elts();
      long *yp = &y.tbl[0][0];
      for (j = 0; j < l; j++) 
         xp[j].LoopHole() = yp[j+lo];
   }
   else {
      for (j = 0; j < l; j++) {
         for (i = 0; i < 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
   

{
   long k, n, i, j, l;
   long NumPrimes = zz_pInfo->NumPrimes;

   long t[4];
   vec_long& s = FFTBuf;

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

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

   long index = zz_pInfo->index;

   if (index >= 0) {
      long *yp = &y.tbl[0][0];
      long q = FFTPrime[index];
      long *Root = &RootTable[index][0];

      FFT(sp, yp, k, q, Root);
      for (j = 0; j < n; j++)
         yp[j] = sp[j];
   }
   else {
      for (i = 0; i < 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);

   if (index >= 0) {
      zz_p *xp = x.elts();
      long *yp = &y.tbl[0][0];
      for (j = 0; j < l; j++) 
         xp[j].LoopHole() = yp[j+lo];
   }
   else {
      for (j = 0; j < l; j++) {
         for (i = 0; i < 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)
{
   long k, n, i, j, l;
   long NumPrimes = zz_pInfo->NumPrimes;

   long t[4];

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

   z.SetSize(k);

   long index = zz_pInfo->index;

   if (index >= 0) {
      long *zp = &z.tbl[0][0];
      long q = FFTPrime[index];
      double qinv = FFTPrimeInv[index];
      long w = TwoInvTable[index][k];
      long *Root = &RootInvTable[index][0];

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

      for (j = 0; j < n; j++) zp[j] = MulMod(zp[j], w, q, qinv);
   }
   else {
      for (i = 0; i < NumPrimes; i++) {
         long *zp = &z.tbl[i][0];
         long q = FFTPrime[i];
         double qinv = FFTPrimeInv[i];
         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);

   if (index >= 0) {
      zz_p *xp = x.rep.elts();
      long *zp = &z.tbl[0][0];
      for (j = 0; j < l; j++) 
         xp[j].LoopHole() = zp[j+lo];
   }
   else {
      for (j = 0; j < l; j++) {
         for (i = 0; i < 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
   

{
   long k, n, i, j;
   long NumPrimes = zz_pInfo->NumPrimes;

   long t[4];
   vec_long& s = FFTBuf;

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

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

   long index = zz_pInfo->index;
   if (index >= 0) {
      long *yp = &y.tbl[0][0];
      long q = FFTPrime[index];
      double qinv = FFTPrimeInv[index];
      long w = TwoInvTable[index][k];
      long *Root = &RootInvTable[index][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 {
            x[j-lo].LoopHole() = y.tbl[0][j];
         }
      }
   }
   else {
      for (i = 0; i < NumPrimes; i++) {
         long *yp = &y.tbl[i][0];
         long q = FFTPrime[i];
         double qinv = FFTPrimeInv[i];
         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)
{
   long k, n, i, j;

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

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

   z.SetSize(k);

   long index = zz_pInfo->index;

   if (index >= 0) {
      long *zp = &z.tbl[0][0];
      const long *xp = &x.tbl[0][0];
      const long *yp = &y.tbl[0][0];
      long q = FFTPrime[index];
      double qinv = FFTPrimeInv[index];

      for (j = 0; j < n; j++)
         zp[j] = MulMod(xp[j], yp[j], q, qinv);
   }
   else {
      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 = FFTPrimeInv[i];
   
         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)
{
   long k, n, i, j;

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

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

   z.SetSize(k);

   long index = zz_pInfo->index;

   if (index >= 0) {
      long *zp = &z.tbl[0][0];
      const long *xp = &x.tbl[0][0];
      const long *yp = &y.tbl[0][0];
      long q = FFTPrime[index];

      for (j = 0; j < n; j++)
         zp[j] = SubMod(xp[j], yp[j], q);
   }
   else {
      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)
{
   long k, n, i, j;

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

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

   z.SetSize(k);

   long index = zz_pInfo->index;

   if (index >= 0) {
      long *zp = &z.tbl[0][0];
      const long *xp = &x.tbl[0][0];
      const long *yp = &y.tbl[0][0];
      long q = FFTPrime[index];

      for (j = 0; j < n; j++)
         zp[j] = AddMod(xp[j], yp[j], q);
   }
   else {
      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
{
   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)
{
   long i, j, l, k, n;

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

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

   long index = zz_pInfo->index;
   
   if (index >= 0) {
      long q = FFTPrime[index];
      const long *ap = &a.tbl[0][0];
      long *xp = &x.tbl[0][0];
      for (j = 0; j < n; j++) {
         long j1 = j << (l-k);
         xp[j1] = AddMod(xp[j1], ap[j], q);
      }
   }
   else {
      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 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_MOD_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_MOD_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_MOD_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: uninitialized modulus");

   if (da <= 2*n-2) {
      rem21(x, a, F);
      return;
   }
   else if (!F.UseFFT || da-n <= NTL_zz_pX_MOD_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("DivRem: uninitialized modulus");

   if (da <= 2*n-2) {
      DivRem21(q, r, a, F);
      return;
   }
   else if (!F.UseFFT || da-n <= NTL_zz_pX_MOD_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;