gf2ex.cpp

数值算法库for Windows

{
   if (a == 0 || b == 0) {
      clear(x);
      return;
   }

   GF2X aa, bb, xx;

   long sx = deg(a) + deg(b) + 1;
   long blocksz = 2*GF2E::degree() - 1;

   if (blocksz >= (1L << (NTL_BITS_PER_LONG-4))/sx)
      Error("overflow in GF2EX KronMul");

   KronSubst(aa, a);
   KronSubst(bb, b);
   mul(xx, aa, bb);

   GF2X c;

   long wc = (blocksz + NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;

   x.rep.SetLength(sx);

   long i;
   for (i = 0; i < sx-1; i++) {
      c.xrep.SetLength(wc);
      ExtractBits(c.xrep.elts(), xx.xrep.elts(), blocksz, i*blocksz);
      c.normalize();
      conv(x.rep[i], c);
   }

   long last_blocksz = deg(xx) - (sx-1)*blocksz + 1;
   wc = (last_blocksz + NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;
   c.xrep.SetLength(wc);

   ExtractBits(c.xrep.elts(), xx.xrep.elts(), last_blocksz, (sx-1)*blocksz);
   c.normalize();
   conv(x.rep[sx-1], c);

   x.normalize();
}



void mul(GF2EX& c, const GF2EX& a, const GF2EX& b)
{
   if (IsZero(a) || IsZero(b)) {
      clear(c);
      return;
   }

   if (&a == &b) {
      sqr(c, a);
      return;
   }

   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sa == 1) {
      mul(c, b, a.rep[0]);
      return;
   }

   if (sb == 1) {
      mul(c, a, b.rep[0]);
      return;
   }

   if (sa < GF2E::KarCross() || sb < GF2E::KarCross()) {
      PlainMul(c, a, b);
      return;
   }

   if (GF2E::WordLength() <= 1) {
      KronMul(c, a, b);
      return;
   }
   

   /* karatsuba */

   long n, hn, sp;

   n = max(sa, sb);
   sp = 0;
   do {
      hn = (n+1) >> 1;
      sp += (hn << 2) - 1;
      n = hn;
   } while (n > 1);

   GF2XVec stk;
   stk.SetSize(sp + 2*(sa+sb)-1, 2*GF2E::WordLength()); 

   long i;

   for (i = 0; i < sa; i++)
      stk[i+sa+sb-1] = rep(a.rep[i]);

   for (i = 0; i < sb; i++)
      stk[i+2*sa+sb-1] = rep(b.rep[i]);

   KarMul(&stk[0], &stk[sa+sb-1], sa, &stk[2*sa+sb-1], sb, 
          &stk[2*(sa+sb)-1]);

   c.rep.SetLength(sa+sb-1);

   for (i = 0; i < sa+sb-1; i++)
      conv(c.rep[i], stk[i]);

   c.normalize();
}


void MulTrunc(GF2EX& x, const GF2EX& a, const GF2EX& b, long n)
{
   GF2EX t;
   mul(t, a, b);
   trunc(x, t, n);
}

void SqrTrunc(GF2EX& x, const GF2EX& a, long n)
{
   GF2EX t;
   sqr(t, a);
   trunc(x, t, n);
}



void PlainDivRem(GF2EX& q, GF2EX& r, 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) {
      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]);
   }

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