lzz_px.cpp

数值算法库for Windows

      for (i = m+1; i >= 1; i--)
         hh[i] = aa[i-1];
      clear(hh[0]);
   }
   else {
      h.rep.SetLength(n);
      hh = h.rep.elts();
      aa = a.rep.elts();
      ff = f.rep.elts();
      negate(z, aa[n-1]);
      if (!IsOne(ff[n]))
         div(z, z, ff[n]);
      for (i = n-1; i >= 1; i--) {
         mul(t, z, ff[i]);
         add(hh[i], aa[i-1], t);
      }
      mul(hh[0], z, ff[0]);
      h.normalize();
   }
}

void MulByXMod(zz_pX& h, const zz_pX& a, const zz_pX& f)
{
   if (&h == &f) {
      zz_pX hh;
      MulByXModAux(hh, a, f);
      h = hh;
   }
   else
      MulByXModAux(h, a, f);
}


void random(zz_pX& x, long n)
{
   long i;

   x.rep.SetLength(n);

   for (i = 0; i < n; i++)
      random(x.rep[i]); 

   x.normalize();
}





void fftRep::SetSize(long NewK)
{
   if (NewK < -1 || NewK >= NTL_BITS_PER_LONG-1)
      Error("bad arg to fftRep::SetSize()");

   if (NewK <= MaxK) {
      k = NewK;
      return;
   }

   if (NumPrimes != zz_pInfo->NumPrimes)
      Error("fftRep: inconsistent use");

   long i, n;
 
   if (MaxK != -1) 
      for (i = 0; i < zz_pInfo->NumPrimes; i++) 
         free(tbl[i]);

   n = 1L << NewK;

   for (i = 0; i < zz_pInfo->NumPrimes; i++) {
      if ( !(tbl[i] = (long *) malloc(n * (sizeof (long)))) )
         Error("out of space in fftRep::SetSize()");
   }

   k = MaxK = NewK;
}

fftRep::fftRep(const fftRep& R)
{
   k = MaxK = R.k;
   NumPrimes = R.NumPrimes;

   if (k < 0) return;

   long i, j, n;

   n = 1L << k;

   for (i = 0; i < NumPrimes; i++) {
      if ( !(tbl[i] = (long *) malloc(n * (sizeof (long)))) )
         Error("out of space in fftRep");

      for (j = 0; j < n; j++)
         tbl[i][j] = R.tbl[i][j];
   }
}

fftRep& fftRep::operator=(const fftRep& R)
{
   if (this == &R) return *this;

   if (NumPrimes != R.NumPrimes)
      Error("fftRep: inconsistent use");

   if (R.k < 0) {
      k = -1;
      return *this;
   }

   if (R.k > MaxK) {
      long i, n;

      if (MaxK != -1) {
         for (i = 0; i < NumPrimes; i++)
            free(tbl[i]);
      }
  
      n = 1L << R.k;
  
      for (i = 0; i < NumPrimes; i++) {
         if ( !(tbl[i] = (long *) malloc(n * (sizeof (long)))) )
            Error("out of space in fftRep");
      }

      k = MaxK = R.k;
   }
   else {
      k = R.k;
   }

   long i, j, n;

   n = 1L << k;

   for (i = 0; i < NumPrimes; i++)
      for (j = 0; j < n; j++)
         tbl[i][j] = R.tbl[i][j];

   return *this;
}



fftRep::~fftRep()
{
   if (MaxK == -1)
      return;

   for (long i = 0; i < NumPrimes; i++)
      free(tbl[i]);
}



static vec_long FFTBuf;




void FromModularRep(zz_p& x, long *a)
{
   long n = zz_pInfo->NumPrimes;
   long p = zz_pInfo->p;
   double pinv = zz_pInfo->pinv;
   long q, s, t;
   long i;
   double y;

#if QUICK_CRT
   y = 0;
   for (i = 0; i < n; i++)
      y = y + ((double) a[i])*zz_pInfo->x[i];

   y = y - long(y*pinv)*p;
   y = y + 0.5;
   while (y >= p) y -= p;
   while (y < 0) y += p;
   q = long(y);
#else
   long Q, r;
   double qq;

   y = 0;
   qq = 0;

   for (i = 0; i < n; i++) {
      r = MulDivRem(Q, a[i], zz_pInfo->u[i], FFTPrime[i], zz_pInfo->x[i]);
      qq = qq + Q;
      y = y + r*FFTPrimeInv[i];
   }

   y = qq + long(y + 0.5);
   y = y - long(y*pinv)*p;
   while (y >= p) y -= p;
   while (y < 0) y += p;
   q = long(y); 

#endif

   t = 0;
   for (i = 0; i < n; i++) {
      s = MulMod(a[i], zz_pInfo->CoeffModP[i], p, pinv);
      t = AddMod(t, s, p);
   }
      

   s = MulMod(q, zz_pInfo->MinusMModP, p, pinv);
   t = AddMod(t, s, p);
   x.LoopHole() = t;
}


void TofftRep(fftRep& y, const zz_pX& x, long k, long lo, long hi)
// computes an n = 2^k point convolution.
// if deg(x) >= 2^k, then x is first reduced modulo X^n-1.
{
   long n, i, j, m, j1;
   vec_long& s = FFTBuf;;
   zz_p accum;
   long NumPrimes = zz_pInfo->NumPrimes;


   if (k > zz_pInfo->MaxRoot) 
      Error("Polynomial too big for FFT");

   if (lo < 0)
      Error("bad arg to TofftRep");

   hi = min(hi, deg(x));

   y.SetSize(k);

   n = 1L << k;

   m = max(hi-lo + 1, 0);

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

   long index = zz_pInfo->index;

   if (index >= 0) {
      for (j = 0; j < n; j++) {
         if (j >= m) {
            y.tbl[0][j] = 0;
         }
         else {
            accum = xx[j+lo];
            for (j1 = j + n; j1 < m; j1 += n)
               add(accum, accum, xx[j1+lo]);
            y.tbl[0][j] = rep(accum);
         }
      }
   }
   else {
      for (j = 0; j < n; j++) {
         if (j >= m) {
            for (i = 0; i < NumPrimes; i++)
               y.tbl[i][j] = 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][j] = t;
            }
         }
      }
   }
   

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

   if (index >= 0) {
      long *Root = &RootTable[index][0];
      long *yp = &y.tbl[0][0];
      FFT(sp, yp, y.k, FFTPrime[index], Root);
      for (j = 0; j < n; j++)
         yp[j] = sp[j];
   } 
   else {
      for (i = 0; i < zz_pInfo->NumPrimes; i++) {
         long *Root = &RootTable[i][0];
         long *yp = &y.tbl[i][0];
         FFT(sp, yp, y.k, FFTPrime[i], Root);
         for (j = 0; j < n; j++)
            yp[j] = sp[j];
      }
   }
}



void RevTofftRep(fftRep& y, const vec_zz_p& x, 
                 long k, long lo, long hi, long offset)
// computes an n = 2^k point convolution of X^offset*x[lo..hi] mod X^n-1
// using "inverted" evaluation points.

{
   long n, i, j, m, j1;
   vec_long& s = FFTBuf;
   zz_p accum;
   long NumPrimes = zz_pInfo->NumPrimes;

   if (k > zz_pInfo->MaxRoot) 
      Error("Polynomial too big for FFT");

   if (lo < 0)
      Error("bad arg to TofftRep");

   hi = min(hi, x.length()-1);

   y.SetSize(k);

   n = 1L << k;

   m = max(hi-lo + 1, 0);

   const zz_p *xx = x.elts();

   long index = zz_pInfo->index;

   offset = offset & (n-1);

   if (index >= 0) {
      for (j = 0; j < n; j++) {
         if (j >= m) {
            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;