zz_px.cpp

数值算法库for Windows

      return;
   }

   ZZ_pX lb;

   if (&q == &b) {
      lb = b;
      bp = lb.rep.elts();
   }
   else
      bp = b.rep.elts();

   if (IsOne(bp[db]))
      LCIsOne = 1;
   else {
      LCIsOne = 0;
      inv(LCInv, bp[db]);
   }

   for (i = 0; i <= da; i++)
      x[i] = rep(a.rep[i]);

   xp = x.elts();

   dq = da - db;
   q.rep.SetLength(dq+1);
   qp = q.rep.elts();

   for (i = dq; i >= 0; i--) {
      conv(t, xp[i+db]);
      if (!LCIsOne)
	 mul(t, t, LCInv);
      qp[i] = t;
      negate(t, t);

      for (j = db-1; j >= 0; j--) {
	 mul(s, rep(t), rep(bp[j]));
	 add(xp[i+j], xp[i+j], s);
      }
   }

   r.rep.SetLength(db);
   for (i = 0; i < db; i++)
      conv(r.rep[i], xp[i]);
   r.normalize();
}


void PlainDiv(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const ZZ_p *bp;
   ZZ_p *qp;
   ZZ *xp;


   ZZ_p LCInv, t;
   static ZZ s;

   da = deg(a);
   db = deg(b);

   if (db < 0) Error("ZZ_pX: division by zero");

   if (da < db) {
      clear(q);
      return;
   }

   ZZ_pX lb;

   if (&q == &b) {
      lb = b;
      bp = lb.rep.elts();
   }
   else
      bp = b.rep.elts();

   if (IsOne(bp[db]))
      LCIsOne = 1;
   else {
      LCIsOne = 0;
      inv(LCInv, bp[db]);
   }

   ZZVec x(da + 1 - db, ZZ_pInfo->ExtendedModulusSize);

   for (i = db; i <= da; i++)
      x[i-db] = rep(a.rep[i]);

   xp = x.elts();

   dq = da - db;
   q.rep.SetLength(dq+1);
   qp = q.rep.elts();

   for (i = dq; i >= 0; i--) {
      conv(t, xp[i]);
      if (!LCIsOne)
	 mul(t, t, LCInv);
      qp[i] = t;
      negate(t, t);

      long lastj = max(0, db-i);

      for (j = db-1; j >= lastj; j--) {
	 mul(s, rep(t), rep(bp[j]));
	 add(xp[i+j-db], xp[i+j-db], s);
      }
   }
}

void PlainRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const ZZ_p *bp;
   ZZ *xp;


   ZZ_p LCInv, t;
   static ZZ s;

   da = deg(a);
   db = deg(b);

   if (db < 0) Error("ZZ_pX: division by zero");

   if (da < db) {
      r = a;
      return;
   }

   bp = b.rep.elts();

   if (IsOne(bp[db]))
      LCIsOne = 1;
   else {
      LCIsOne = 0;
      inv(LCInv, bp[db]);
   }

   ZZVec x(da + 1, ZZ_pInfo->ExtendedModulusSize);

   for (i = 0; i <= da; i++)
      x[i] = rep(a.rep[i]);

   xp = x.elts();

   dq = da - db;

   for (i = dq; i >= 0; i--) {
      conv(t, xp[i+db]);
      if (!LCIsOne)
	 mul(t, t, LCInv);
      negate(t, t);

      for (j = db-1; j >= 0; j--) {
	 mul(s, rep(t), rep(bp[j]));
	 add(xp[i+j], xp[i+j], s);
      }
   }

   r.rep.SetLength(db);
   for (i = 0; i < db; i++)
      conv(r.rep[i], xp[i]);
   r.normalize();
}

void mul(ZZ_pX& x, const ZZ_pX& a, const ZZ_p& b)
{
   if (IsZero(b)) {
      clear(x);
      return;
   }

   if (IsOne(b)) {
      x = a;
      return;
   }

   ZZ_pTemp TT; ZZ_p& t = TT.val();

   long i, da;

   const ZZ_p *ap;
   ZZ_p* xp;


   t = b;

   da = deg(a);
   x.rep.SetLength(da+1);
   ap = a.rep.elts();
   xp = x.rep.elts();

   for (i = 0; i <= da; i++) 
      mul(xp[i], ap[i], t);

   x.normalize();
}

void mul(ZZ_pX& x, const ZZ_pX& a, long b)
{
   ZZ_pTemp TT;  ZZ_p& T = TT.val();
   conv(T, b);
   mul(x, a, T);
}


void PlainGCD(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b)
{
   ZZ_p t;

   if (IsZero(b))
      x = a;
   else if (IsZero(a))
      x = b;
   else {
      long n = max(deg(a),deg(b)) + 1;
      ZZ_pX u(INIT_SIZE, n), v(INIT_SIZE, n);
      ZZVec tmp(n, ZZ_pInfo->ExtendedModulusSize);

      u = a;
      v = b;
      do {
         PlainRem(u, u, v, tmp);
         swap(u, v);
      } while (!IsZero(v));

      x = u;
   }

   if (IsZero(x)) return;
   if (IsOne(LeadCoeff(x))) return;

   /* make gcd monic */


   inv(t, LeadCoeff(x)); 
   mul(x, x, t); 
}



         

void PlainXGCD(ZZ_pX& d, ZZ_pX& s, ZZ_pX& t, const ZZ_pX& a, const ZZ_pX& b)
{
   ZZ_p z;


   if (IsZero(b)) {
      set(s);
      clear(t);
      d = a;
   }
   else if (IsZero(a)) {
      clear(s);
      set(t);
      d = b;
   }
   else {
      long e = max(deg(a), deg(b)) + 1;

      ZZ_pX temp(INIT_SIZE, e), u(INIT_SIZE, e), v(INIT_SIZE, e), 
            u0(INIT_SIZE, e), v0(INIT_SIZE, e), 
            u1(INIT_SIZE, e), v1(INIT_SIZE, e), 
            u2(INIT_SIZE, e), v2(INIT_SIZE, e), q(INIT_SIZE, e);


      set(u1); clear(v1);
      clear(u2); set(v2);
      u = a; v = b;

      do {
         DivRem(q, u, u, v);
         swap(u, v);
         u0 = u2;
         v0 = v2;
         mul(temp, q, u2);
         sub(u2, u1, temp);
         mul(temp, q, v2);
         sub(v2, v1, temp);
         u1 = u0;
         v1 = v0;
      } while (!IsZero(v));

      d = u;
      s = u1;
      t = v1;
   }

   if (IsZero(d)) return;
   if (IsOne(LeadCoeff(d))) return;

   /* make gcd monic */

   inv(z, LeadCoeff(d));
   mul(d, d, z);
   mul(s, s, z);
   mul(t, t, z);
}


void MulMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, const ZZ_pX& f)
{
   if (deg(a) >= deg(f) || deg(b) >= deg(f) || deg(f) == 0) 
      Error("MulMod: bad args");

   ZZ_pX t;

   mul(t, a, b);
   rem(x, t, f);
}

void SqrMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("SqrMod: bad args");

   ZZ_pX t;

   sqr(t, a);
   rem(x, t, f);
}


void InvMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args");

   ZZ_pX d, t;

   XGCD(d, x, t, a, f);
   if (!IsOne(d))
      Error("ZZ_pX InvMod: can't compute multiplicative inverse");
}

long InvModStatus(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvModStatus: bad args");
   ZZ_pX d, t;

   XGCD(d, x, t, a, f);
   if (!IsOne(d)) {
      x = d;
      return 1;
   }
   else
      return 0;
}




static
void MulByXModAux(ZZ_pX& h, const ZZ_pX& a, const ZZ_pX& f)
{
   long i, n, m;
   ZZ_p* hh;
   const ZZ_p *aa, *ff;

   ZZ_p t, z;

   n = deg(f);
   m = deg(a);

   if (m >= n || n == 0) Error("MulByXMod: bad args");

   if (m < 0) {
      clear(h);
      return;
   }

   if (m < n-1) {
      h.rep.SetLength(m+2);
      hh = h.rep.elts();
      aa = a.rep.elts();
      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;
   }

   ZZ_pInfo->check();

   if (MaxK == -1)
      NumPrimes = ZZ_pInfo->NumPrimes;
   else {
      if (NumPrimes != ZZ_pInfo->NumPrimes)
         Error("FFTRep: inconsistent use");
   }

   long i, n;

   if (MaxK == -1) {
      tbl = (long **) malloc(NumPrimes * (sizeof (long *)) );
      if (!tbl)
         Error("out of space in FFTRep::SetSize()");
   }
   else {
      for (i = 0; i < NumPrimes; i++) 
         free(tbl[i]);
   }

   n = 1L << NewK;

   for (i = 0; i < 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;
   tbl = 0;
   NumPrimes = 0;

   if (k < 0) return;

   NumPrimes = R.NumPrimes;

   long i, j, n;
 
   tbl = (long **) malloc(NumPrimes * (sizeof (long *)) );
   if (!tbl)
      Error("out of space in FFTRep");

   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 (MaxK >= 0 && R.MaxK >= 0 && NumPrimes != R.NumPrimes)
      Error("FFTRep: inconsistent use");

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

   NumPrimes = R.NumPrimes;

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

      if (MaxK == -1) {
         tbl = (long **) malloc(NumPrimes * (sizeof (long *)) );
         if (!tbl)
            Error("out of space in FFTRep");
      }
      else {
         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]);

   free(tbl);
}



void ZZ_pXModRep::SetSize(long NewN)
{
   ZZ_pInfo->check();

   NumPrimes = ZZ_pInfo->NumPrimes;

   if (NewN < 0)
      Error("bad arg to ZZ_pXModRep::SetSize()");

   if (NewN <= MaxN) {
      n = NewN;
      return;
   }

   long i;
 

   if (MaxN == 0) {
      tbl = (long **) malloc(ZZ_pInfo->NumPrimes * (sizeof (long *)) );
      if (!tbl)
         Error("out of space in ZZ_pXModRep::SetSize()");
   }
   else {
      for (i = 0; i < ZZ_pInfo->NumPrimes; i++) 
         free(tbl[i]);
   }

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

   n = MaxN = NewN;
}

ZZ_pXModRep::~ZZ_pXModRep()
{
   if (MaxN == 0)
      return;

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

   free(tbl);
}


static vec_long ModularRepBuf;
static vec_long FFTBuf;




void ToModularRep(vec_long& x, const ZZ_p& a)
{
   ZZ_pInfo->check();
   ZZ_p_rem_struct_eval(ZZ_pInfo->rem_struct, &x[0], rep(a));
}


// NOTE: earlier versions used Kahan summation...
// we no longer do this, as it is less portable than I thought.



void FromModularRep(ZZ_p& x, const vec_long& a)
{
   ZZ_pInfo->check();

   long n = ZZ_pInfo->NumPrimes;
   static ZZ q, s, t;
   long i;
   double y;

   if (ZZ_p_crt_struct_special(ZZ_pInfo->crt_struct)) {
      ZZ_p_crt_struct_eval(ZZ_pInfo->crt_struct, t, &a[0]);
      x.LoopHole() = t;
      return;
   }
      

   if (ZZ_pInfo->QuickCRT) {
      y = 0;
      for (i = 0; i < n; i++)
         y += ((double) a[i])*ZZ_pInfo->x[i];

      conv(q, (y + 0.5)); 
   } else {
      long Q, r;
      static ZZ qq;

      y = 0;

      clear(q);

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

      conv(qq, (y + 0.5));
      add(q, q, qq);
   }

   ZZ_p_crt_struct_eval(ZZ_pInfo->crt_struct, t, &a[0]);

   mul(s, q, ZZ_pInfo->MinusMModP);
   add(t, t, s);

   conv(x, 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.
{
   ZZ_pInfo->check();

   long n, i, j, m, j1;
   vec_long& t = ModularRepBuf;
   vec_long& s = FFTBuf;;
   ZZ_p accum;


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

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

   t.SetLength(ZZ_pInfo->NumPrimes);

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

   y.SetSize(k);

   n = 1L << k;

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

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

   for (j = 0; j < n; j++) {
      if (j >= m) {
         for (i = 0; i < ZZ_pInfo->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]);
         ToModularRep(t, accum);
         for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
            y.tbl[i][j] = t[i];
         }
      }
   }


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

   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.

{
   ZZ_pInfo->check();

   long n, i, j, m, j1;
   vec_long& t = ModularRepBuf;
   vec_long& s = FFTBuf;
   ZZ_p accum;

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

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

   t.SetLength(ZZ_pInfo->NumPrimes);

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

   y.SetSize(k);

   n = 1L << k;

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

   const ZZ_p *xx = x.elts();

   offset = offset & (n-1);

   for (j = 0; j < n; j++) {
      if (j >= m) {
         for (i = 0; i < ZZ_pInfo->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]);
         ToModularRep(t, accum);
         for (i = 0; i < ZZ_pInfo->NumPrimes; i++) {
            y.tbl[i][offset] = t[i];

         }
      }

      offset = (offset + 1) & (n-1);
   }


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

   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);
   }