lzz_pex.c

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

{
   long n = a.rep.length();
   x.rep.SetLength(n);

   const zz_pE* ap = a.rep.elts();
   zz_pE* xp = x.rep.elts();
   long i;

   for (i = n; i; i--, ap++, xp++)
      negate((*xp), (*ap));
}



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

   zz_pE 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_pEX& h, const zz_pEX& a, const zz_pEX& f)
{
   if (&h == &f) {
      zz_pEX hh;
      MulByXModAux(hh, a, f);
      h = hh;
   }
   else
      MulByXModAux(h, a, f);
}



void PlainMul(zz_pEX& x, const zz_pEX& a, const zz_pEX& b)
{
   long da = deg(a);
   long db = deg(b);

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

   long d = da+db;



   const zz_pE *ap, *bp;
   zz_pE *xp;
   
   zz_pEX la, lb;

   if (&x == &a) {
      la = a;
      ap = la.rep.elts();
   }
   else
      ap = a.rep.elts();

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

   x.rep.SetLength(d+1);

   xp = x.rep.elts();

   long i, j, jmin, jmax;
   static zz_pX t, accum;

   for (i = 0; i <= d; i++) {
      jmin = max(0, i-db);
      jmax = min(da, i);
      clear(accum);
      for (j = jmin; j <= jmax; j++) {
	 mul(t, rep(ap[j]), rep(bp[i-j]));
	 add(accum, accum, t);
      }
      conv(xp[i], accum);
   }
   x.normalize();
}

void SetSize(vec_zz_pX& x, long n, long m)
{
   x.SetLength(n);
   long i;
   for (i = 0; i < n; i++)
      x[i].rep.SetMaxLength(m);
}



void PlainDivRem(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const zz_pE *bp;
   zz_pE *qp;
   zz_pX *xp;


   zz_pE LCInv, t;
   zz_pX s;

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

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

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

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

   vec_zz_pX x;

   SetSize(x, da+1, 2*zz_pE::degree());

   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 PlainRem(zz_pEX& r, const zz_pEX& a, const zz_pEX& b, vec_zz_pX& x)
{
   long da, db, dq, i, j, LCIsOne;
   const zz_pE *bp;
   zz_pX *xp;


   zz_pE LCInv, t;
   zz_pX s;

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

   if (db < 0) Error("zz_pEX: 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]);
   }

   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 PlainDivRem(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEX& b, 
     vec_zz_pX& x)
{
   long da, db, dq, i, j, LCIsOne;
   const zz_pE *bp;
   zz_pE *qp;
   zz_pX *xp;


   zz_pE LCInv, t;
   zz_pX s;

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

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

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

   zz_pEX 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_pEX& q, const zz_pEX& a, const zz_pEX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const zz_pE *bp;
   zz_pE *qp;
   zz_pX *xp;


   zz_pE LCInv, t;
   zz_pX s;

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

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

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

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

   vec_zz_pX x;
   SetSize(x, da+1-db, 2*zz_pE::degree());

   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_pEX& r, const zz_pEX& a, const zz_pEX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const zz_pE *bp;
   zz_pX *xp;


   zz_pE LCInv, t;
   zz_pX s;

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

   if (db < 0) Error("zz_pEX: 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]);
   }

   vec_zz_pX x;
   SetSize(x, da + 1, 2*zz_pE::degree());

   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 RightShift(zz_pEX& x, const zz_pEX& a, long n)
{
   if (n < 0) {
      if (n < -NTL_MAX_LONG) Error("overflow in RightShift");
      LeftShift(x, a, -n);
      return;
   }

   long da = deg(a);
   long i;
 
   if (da < n) {
      clear(x);
      return;
   }

   if (&x != &a)
      x.rep.SetLength(da-n+1);

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

   if (&x == &a)
      x.rep.SetLength(da-n+1);

   x.normalize();
}

void LeftShift(zz_pEX& x, const zz_pEX& a, long n)
{
   if (n < 0) {
      if (n < -NTL_MAX_LONG) Error("overflow in LeftShift");
      RightShift(x, a, -n);
      return;
   }

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

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

   long m = a.rep.length();

   x.rep.SetLength(m+n);

   long i;
   for (i = m-1; i >= 0; i--)
      x.rep[i+n] = a.rep[i];

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



void NewtonInv(zz_pEX& c, const zz_pEX& a, long e)
{
   zz_pE x;

   inv(x, ConstTerm(a));

   if (e == 1) {
      conv(c, x);
      return;
   }

   static vec_long E;
   E.SetLength(0);
   append(E, e);
   while (e > 1) {
      e = (e+1)/2;
      append(E, e);
   }

   long L = E.length();

   zz_pEX g, g0, g1, g2;


   g.rep.SetMaxLength(e);
   g0.rep.SetMaxLength(e);
   g1.rep.SetMaxLength((3*e+1)/2);
   g2.rep.SetMaxLength(e);

   conv(g, x);

   long i;

   for (i = L-1; i > 0; i--) {
      // lift from E[i] to E[i-1]

      long k = E[i];
      long l = E[i-1]-E[i];

      trunc(g0, a, k+l);

      mul(g1, g0, g);
      RightShift(g1, g1, k);
      trunc(g1, g1, l);

      mul(g2, g1, g);
      trunc(g2, g2, l);
      LeftShift(g2, g2, k);

      sub(g, g, g2);
   }

   c = g;
}

void InvTrunc(zz_pEX& c, const zz_pEX& a, long e)
{
   if (e < 0) Error("InvTrunc: bad args");
   if (e == 0) {
      clear(c);
      return;
   }

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

   NewtonInv(c, a, e);
}




const long zz_pEX_MOD_PLAIN = 0;
const long zz_pEX_MOD_MUL = 1;


void build(zz_pEXModulus& F, const zz_pEX& f)
{
   long n = deg(f);

   if (n <= 0) Error("build(zz_pEXModulus,zz_pEX): deg(f) <= 0");

   if (n >= (1L << (NTL_BITS_PER_LONG-4))/zz_pE::degree())
      Error("build(zz_pEXModulus,zz_pEX): overflow");

   F.tracevec.SetLength(0);

   F.f = f;
   F.n = n;

   if (F.n < zz_pE::ModCross()) {
      F.method = zz_pEX_MOD_PLAIN;
   }
   else {
      F.method = zz_pEX_MOD_MUL;
      zz_pEX P1;
      zz_pEX P2;

      CopyReverse(P1, f, n);
      InvTrunc(P2, P1, n-1);
      CopyReverse(P1, P2, n-2);
      trunc(F.h0, P1, n-2);
      trunc(F.f0, f, n);
      F.hlc = ConstTerm(P2);
   }
}



zz_pEXModulus::zz_pEXModulus()
{
   n = -1;
   method = zz_pEX_MOD_PLAIN;
}


zz_pEXModulus::~zz_pEXModulus() 
{ 
}



zz_pEXModulus::zz_pEXModulus(const zz_pEX& ff)
{
   n = -1;
   method = zz_pEX_MOD_PLAIN;

   build(*this, ff);
}


void UseMulRem21(zz_pEX& r, const zz_pEX& a, const zz_pEXModulus& F)
{
   zz_pEX P1;
   zz_pEX P2;

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
   add(P2, P2, P1);
   mul(P1, P2, F.f0);
   trunc(P1, P1, F.n);
   trunc(r, a, F.n);
   sub(r, r, P1);
}

void UseMulDivRem21(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEXModulus& F)
{
   zz_pEX P1;
   zz_pEX P2;

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
   add(P2, P2, P1);
   mul(P1, P2, F.f0);
   trunc(P1, P1, F.n);
   trunc(r, a, F.n);
   sub(r, r, P1);
   q = P2;
}

void UseMulDiv21(zz_pEX& q, const zz_pEX& a, const zz_pEXModulus& F)
{
   zz_pEX P1;
   zz_pEX P2;

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
   add(P2, P2, P1);
   q = P2;

}


void rem(zz_pEX& x, const zz_pEX& a, const zz_pEXModulus& F)
{
   if (F.method == zz_pEX_MOD_PLAIN) {
      PlainRem(x, a, F.f);
      return;
   }

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

   if (da <= 2*n-2) {
      UseMulRem21(x, a, F);
      return;
   }

   zz_pEX 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();

      UseMulRem21(buf, buf, F);

      a_len -= amt;
   }

   x = buf;
}

void DivRem(zz_pEX& q, zz_pEX& r, const zz_pEX& a, const zz_pEXModulus& F)
{
   if (F.method == zz_pEX_MOD_PLAIN) {
      PlainDivRem(q, r, a, F.f);
      return;
   }

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

   if (da <= 2*n-2) {
      UseMulDivRem21(q, r, a, F);
      return;
   }

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

   zz_pEX 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();

      UseMulDivRem21(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_pEX& q, const zz_pEX& a, const zz_pEXModulus& F)
{
   if (F.method == zz_pEX_MOD_PLAIN) {
      PlainDiv(q, a, F.f);
      return;
   }

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

   if (da <= 2*n-2) {
      UseMulDiv21(q, a, F);
      return;
   }

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

   zz_pEX 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)
         UseMulDivRem21(qbuf, buf, buf, F);
      else
         UseMulDiv21(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_pEX& c, const zz_pEX& a, const zz_pEX& b, const zz_pEXModulus& F)
{