gf2ex.c

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

   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;

      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(GF2EX& r, const GF2EX& a, const GF2EX& b, GF2XVec& x)
{
   long da, db, dq, i, j, LCIsOne;
   const GF2E *bp;
   GF2X *xp;


   GF2E LCInv, t;
   GF2X s;

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

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

      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(GF2EX& q, GF2EX& r, const GF2EX& a, const GF2EX& b, GF2XVec& x)
{
   long da, db, dq, i, j, LCIsOne;
   const GF2E *bp;
   GF2E *qp;
   GF2X *xp;


   GF2E LCInv, t;
   GF2X s;

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

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

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

   GF2EX 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;

      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(GF2EX& q, const GF2EX& a, const GF2EX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const GF2E *bp;
   GF2E *qp;
   GF2X *xp;


   GF2E LCInv, t;
   GF2X s;

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

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

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

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

   GF2XVec x(da + 1 - db, 2*GF2E::WordLength());

   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;

      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(GF2EX& r, const GF2EX& a, const GF2EX& b)
{
   long da, db, dq, i, j, LCIsOne;
   const GF2E *bp;
   GF2X *xp;


   GF2E LCInv, t;
   GF2X s;

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

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

   GF2XVec x(da + 1, 2*GF2E::WordLength());

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

      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(GF2EX& x, const GF2EX& a, const GF2E& b)
{
   if (IsZero(a) || IsZero(b)) {
      clear(x);
      return;
   }

   GF2X bb, t;
   long i, da;

   const GF2E *ap;
   GF2E* xp;

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

   for (i = 0; i <= da; i++) {
      mul(t, rep(ap[i]), bb);
      conv(xp[i], t);
   }

   x.normalize();
}

void mul(GF2EX& x, const GF2EX& a, GF2 b)
{
   if (b == 0)
      clear(x);
   else
      x = a;
}

void mul(GF2EX& x, const GF2EX& a, long b)
{
   if ((b & 1) == 0)
      clear(x);
   else
      x = a;
}


void GCD(GF2EX& x, const GF2EX& a, const GF2EX& b)
{
   GF2E t;

   if (IsZero(b))
      x = a;
   else if (IsZero(a))
      x = b;
   else {
      long n = max(deg(a),deg(b)) + 1;
      GF2EX u(INIT_SIZE, n), v(INIT_SIZE, n);
      GF2XVec tmp(n, 2*GF2E::WordLength());

      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 XGCD(GF2EX& d, GF2EX& s, GF2EX& t, const GF2EX& a, const GF2EX& b)
{
   GF2E 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;

      GF2EX 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);
         add(u2, u1, temp);
         mul(temp, q, v2);
         add(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(GF2EX& x, const GF2EX& a, const GF2EX& b, const GF2EX& f)
{
   if (deg(a) >= deg(f) || deg(b) >= deg(f) || deg(f) == 0) 
      Error("MulMod: bad args");

   GF2EX t;

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

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

   GF2EX t;

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


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

   GF2EX d, t;

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

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

   GF2EX d, t;

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




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

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




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

   x.rep.SetLength(n);

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

   x.normalize();
}


void CopyReverse(GF2EX& x, const GF2EX& a, long hi)

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

{
   long i, j, n, m;

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

   x.rep.SetLength(n);

   const GF2E* ap = a.rep.elts();
   GF2E* 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 trunc(GF2EX& x, const GF2EX& a, long m)

// x = a % X^m, output may alias input 

{
   if (m < 0) Error("trunc: bad args");

   if (&x == &a) {
      if (x.rep.length() > m) {
         x.rep.SetLength(m);
         x.normalize();
      }
   }
   else {
      long n;
      long i;
      GF2E* xp;
      const GF2E* ap;

      n = min(a.rep.length(), m);
      x.rep.SetLength(n);

      xp = x.rep.elts();
      ap = a.rep.elts();

      for (i = 0; i < n; i++) xp[i] = ap[i];

      x.normalize();
   }
}

void NewtonInvTrunc(GF2EX& c, const GF2EX& a, long e)
{
   GF2E 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();

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

      add(g, g, g2);
   }

   c = g;
}


void InvTrunc(GF2EX& c, const GF2EX& 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");

   NewtonInvTrunc(c, a, e);
}



const long GF2EX_MOD_PLAIN = 0;
const long GF2EX_MOD_MUL = 1;

void build(GF2EXModulus& F, const GF2EX& f)
{
   long n = deg(f);

   if (n <= 0) Error("build(GF2EXModulus,GF2EX): deg(f) <= 0");

   if (n >= (1L << (NTL_BITS_PER_LONG-4))/GF2E::degree())
      Error("build(GF2EXModulus,GF2EX): overflow");

   F.tracevec.SetLength(0);

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

   if (F.n < GF2E::ModCross()) {
      F.method = GF2EX_MOD_PLAIN;
   }
   else {
      F.method = GF2EX_MOD_MUL;
      GF2EX P1;
      GF2EX 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);
   }
}

GF2EXModulus::GF2EXModulus()
{
   n = -1;
   method = GF2EX_MOD_PLAIN;
}



GF2EXModulus::GF2EXModulus(const GF2EX& ff)
{
   n = -1;
   method = GF2EX_MOD_PLAIN;

   build(*this, ff);
}





void UseMulRem21(GF2EX& r, const GF2EX& a, const GF2EXModulus& F)
{
   GF2EX P1;
   GF2EX 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);
   add(r, r, P1);
}

void UseMulDivRem21(GF2EX& q, GF2EX& r, const GF2EX& a, const GF2EXModulus& F)
{
   GF2EX P1;
   GF2EX 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);
   add(r, r, P1);
   q = P2;
}

void UseMulDiv21(GF2EX& q, const GF2EX& a, const GF2EXModulus& F)
{
   GF2EX P1;
   GF2EX 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(GF2EX& x, const GF2EX& a, const GF2EXModulus& F)
{
   if (F.method == GF2EX_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;
   }

   GF2EX 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(GF2EX& q, GF2EX& r, const GF2EX& a, const GF2EXModulus& F)
{
   if (F.method == GF2EX_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;
   }

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

   GF2EX 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(GF2EX& q, const GF2EX& a, const GF2EXModulus& F)
{
   if (F.method == GF2EX_MOD_PLAIN) {