gf2ex.c

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

#include <NTL/GF2EX.h>
#include <NTL/vec_vec_GF2.h>

#include <NTL/new.h>

NTL_START_IMPL



const GF2EX& GF2EX::zero()
{
   static GF2EX z;
   return z;
}



istream& operator>>(istream& s, GF2EX& x)
{
   s >> x.rep;
   x.normalize();
   return s;
}

ostream& operator<<(ostream& s, const GF2EX& a)
{
   return s << a.rep;
}


void GF2EX::normalize()
{
   long n;
   const GF2E* p;

   n = rep.length();
   if (n == 0) return;
   p = rep.elts() + (n-1);
   while (n > 0 && IsZero(*p)) {
      p--; 
      n--;
   }
   rep.SetLength(n);
}


long IsZero(const GF2EX& a)
{
   return a.rep.length() == 0;
}


long IsOne(const GF2EX& a)
{
    return a.rep.length() == 1 && IsOne(a.rep[0]);
}

void GetCoeff(GF2E& x, const GF2EX& a, long i)
{
   if (i < 0 || i > deg(a))
      clear(x);
   else
      x = a.rep[i];
}

void SetCoeff(GF2EX& x, long i, const GF2E& a)
{
   long j, m;

   if (i < 0) 
      Error("SetCoeff: negative index");

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

   m = deg(x);

   if (i > m) {
      long pos = x.rep.position(a);
      x.rep.SetLength(i+1);

      if (pos != -1)
         x.rep[i] = x.rep.RawGet(pos);
      else
         x.rep[i] = a;

      for (j = m+1; j < i; j++)
         clear(x.rep[j]);
   }
   else
      x.rep[i] = a;

   x.normalize();
}

void SetCoeff(GF2EX& x, long i, GF2 a)
{
   if (i < 0)
      Error("SetCoeff: negative index");

   if (a == 1)
      SetCoeff(x, i);
   else
      SetCoeff(x, i, GF2E::zero());
}

void SetCoeff(GF2EX& x, long i, long a)
{
   if (i < 0)
      Error("SetCoeff: negative index");

   if ((a & 1) == 1)
      SetCoeff(x, i);
   else
      SetCoeff(x, i, GF2E::zero());
}

void SetCoeff(GF2EX& x, long i)
{
   long j, m;

   if (i < 0) 
      Error("coefficient index out of range");

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

   m = deg(x);

   if (i > m) {
      x.rep.SetLength(i+1);
      for (j = m+1; j < i; j++)
         clear(x.rep[j]);
   }
   set(x.rep[i]);
   x.normalize();
}


void SetX(GF2EX& x)
{
   clear(x);
   SetCoeff(x, 1);
}


long IsX(const GF2EX& a)
{
   return deg(a) == 1 && IsOne(LeadCoeff(a)) && IsZero(ConstTerm(a));
}
      
      

const GF2E& coeff(const GF2EX& a, long i)
{
   if (i < 0 || i > deg(a))
      return GF2E::zero();
   else
      return a.rep[i];
}


const GF2E& LeadCoeff(const GF2EX& a)
{
   if (IsZero(a))
      return GF2E::zero();
   else
      return a.rep[deg(a)];
}

const GF2E& ConstTerm(const GF2EX& a)
{
   if (IsZero(a))
      return GF2E::zero();
   else
      return a.rep[0];
}



void conv(GF2EX& x, const GF2E& a)
{
   if (IsZero(a))
      x.rep.SetLength(0);
   else {
      x.rep.SetLength(1);
      x.rep[0] = a;
   }
}

void conv(GF2EX& x, long a)
{
   if (a & 1)
      set(x);
   else
      clear(x);
}

void conv(GF2EX& x, GF2 a)
{
   if (a == 1)
      set(x);
   else
      clear(x);
}

void conv(GF2EX& x, const ZZ& a)
{
   if (IsOdd(a))
      set(x);
   else
      clear(x);
}

void conv(GF2EX& x, const GF2X& aa)
{
   GF2X a = aa; // in case a aliases the rep of a coefficient of x
   
   long n = deg(a)+1;
   long i;

   x.rep.SetLength(n);
   for (i = 0; i < n; i++)
      conv(x.rep[i], coeff(a, i));
}

void conv(GF2EX& x, const vec_GF2E& a)
{
   x.rep = a;
   x.normalize();
}


void add(GF2EX& x, const GF2EX& a, const GF2EX& b)
{
   long da = deg(a);
   long db = deg(b);
   long minab = min(da, db);
   long maxab = max(da, db);
   x.rep.SetLength(maxab+1);

   long i;
   const GF2E *ap, *bp; 
   GF2E* xp;

   for (i = minab+1, ap = a.rep.elts(), bp = b.rep.elts(), xp = x.rep.elts();
        i; i--, ap++, bp++, xp++)
      add(*xp, (*ap), (*bp));

   if (da > minab && &x != &a)
      for (i = da-minab; i; i--, xp++, ap++)
         *xp = *ap;
   else if (db > minab && &x != &b)
      for (i = db-minab; i; i--, xp++, bp++)
         *xp = *bp;
   else
      x.normalize();
}

void add(GF2EX& x, const GF2EX& a, const GF2E& b)
{
   long n = a.rep.length();
   if (n == 0) {
      conv(x, b);
   }
   else if (&x == &a) {
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else if (x.rep.MaxLength() == 0) {
      x = a;
      add(x.rep[0], a.rep[0], b);
      x.normalize();
   }
   else {
      // ugly...b could alias a coeff of x

      GF2E *xp = x.rep.elts();
      add(xp[0], a.rep[0], b);
      x.rep.SetLength(n);
      xp = x.rep.elts();
      const GF2E *ap = a.rep.elts();
      long i;
      for (i = 1; i < n; i++)
         xp[i] = ap[i];
      x.normalize();
   }
}

void add(GF2EX& x, const GF2EX& a, GF2 b)
{
   if (a.rep.length() == 0) {
      conv(x, b);
   }
   else {
      if (&x != &a) x = a;
      add(x.rep[0], x.rep[0], b);
      x.normalize();
   }
}

void add(GF2EX& x, const GF2EX& a, long b)
{
   if (a.rep.length() == 0) {
      conv(x, b);
   }
   else {
      if (&x != &a) x = a;
      add(x.rep[0], x.rep[0], b);
      x.normalize();
   }
}


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

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

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

   long d = da+db;

   const GF2E *ap, *bp;
   GF2E *xp;
   
   GF2EX 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;
   GF2X 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 sqr(GF2EX& x, const GF2EX& a)
{
   long da = deg(a);

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

   x.rep.SetLength(2*(da+1));
   long i;

   for (i = da; i > 0; i--) {
      sqr(x.rep[2*i], a.rep[i]);
      clear(x.rep[2*i-1]);
   }

   sqr(x.rep[0], a.rep[0]);

   x.normalize();
}


#if 0

static
void PlainMul(GF2X *xp, const GF2X *ap, long sa, const GF2X *bp, long sb)
{
   if (sa == 0 || sb == 0) return;

   long sx = sa+sb-1;

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

   for (i = 0; i < sx; i++) {
      jmin = max(0, i-sb+1);
      jmax = min(sa-1, i);
      clear(accum);
      for (j = jmin; j <= jmax; j++) {
         mul(t, ap[j], bp[i-j]);
         add(accum, accum, t);
      }
      xp[i] = accum;
   }
}

#endif

static void PlainMul1(GF2X *xp, const GF2X *ap, long sa, const GF2X& b)
{
   long i;

   for (i = 0; i < sa; i++)
      mul(xp[i], ap[i], b);
}




inline
void q_add(GF2X& x, const GF2X& a, const GF2X& b)

// This is a quick-and-dirty add rotine used by the karatsuba routine.
// It assumes that the output already has enough space allocated,
// thus avoiding any procedure calls.
// WARNING: it also accesses the underlying WordVector representation
// directly...that is dirty!.
// It shaves a few percent off the running time.

{
   _ntl_ulong *xp = x.xrep.elts();
   const _ntl_ulong *ap = a.xrep.elts();
   const _ntl_ulong *bp = b.xrep.elts();

   long sa = ap[-1];
   long sb = bp[-1];

   long i;

   if (sa == sb) {
      for (i = 0; i < sa; i++)
         xp[i] = ap[i] ^ bp[i];

      i = sa-1;
      while (i >= 0 && !xp[i]) i--;
      xp[-1] = i+1;
   }
   else if (sa < sb) {
      for (i = 0; i < sa; i++)
         xp[i] = ap[i] ^ bp[i];

      for (; i < sb; i++)
         xp[i] = bp[i];

      xp[-1] = sb;
   }
   else { // sa > sb
      for (i = 0; i < sb; i++)
         xp[i] = ap[i] ^ bp[i];

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

      xp[-1] = sa;
   }
}


inline
void q_copy(GF2X& x, const GF2X& a)
// see comments for q_add above

{
   _ntl_ulong *xp = x.xrep.elts();
   const _ntl_ulong *ap = a.xrep.elts();

   long sa = ap[-1];
   long i;

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

   xp[-1] = sa;
}



static
void KarFold(GF2X *T, const GF2X *b, long sb, long hsa)
{
   long m = sb - hsa;
   long i;

   for (i = 0; i < m; i++)
      q_add(T[i], b[i], b[hsa+i]);

   for (i = m; i < hsa; i++)
      q_copy(T[i], b[i]);
}


static
void KarAdd(GF2X *T, const GF2X *b, long sb)
{
   long i;

   for (i = 0; i < sb; i++)
      q_add(T[i], T[i], b[i]);
}

static
void KarFix(GF2X *c, const GF2X *b, long sb, long hsa)
{
   long i;

   for (i = 0; i < hsa; i++)
      q_copy(c[i], b[i]);

   for (i = hsa; i < sb; i++)
      q_add(c[i], c[i], b[i]);
}



static
void KarMul(GF2X *c, const GF2X *a, 
            long sa, const GF2X *b, long sb, GF2X *stk)
{
   if (sa < sb) {
      { long t = sa; sa = sb; sb = t; }
      { const GF2X *t = a; a = b; b = t; }
   }

   if (sb == 1) {  
      if (sa == 1) 
         mul(*c, *a, *b);
      else
         PlainMul1(c, a, sa, *b);

      return;
   }

   if (sb == 2 && sa == 2) {
      mul(c[0], a[0], b[0]);
      mul(c[2], a[1], b[1]);
      q_add(stk[0], a[0], a[1]);
      q_add(stk[1], b[0], b[1]);
      mul(c[1], stk[0], stk[1]);
      q_add(c[1], c[1], c[0]);
      q_add(c[1], c[1], c[2]);
      
      return;
   }

   long hsa = (sa + 1) >> 1;

   if (hsa < sb) {
      /* normal case */

      long hsa2 = hsa << 1;

      GF2X *T1, *T2, *T3;

      T1 = stk; stk += hsa;
      T2 = stk; stk += hsa;
      T3 = stk; stk += hsa2 - 1;

      /* compute T1 = a_lo + a_hi */

      KarFold(T1, a, sa, hsa);

      /* compute T2 = b_lo + b_hi */

      KarFold(T2, b, sb, hsa);

      /* recursively compute T3 = T1 * T2 */

      KarMul(T3, T1, hsa, T2, hsa, stk);

      /* recursively compute a_hi * b_hi into high part of c */
      /* and subtract from T3 */

      KarMul(c + hsa2, a+hsa, sa-hsa, b+hsa, sb-hsa, stk);
      KarAdd(T3, c + hsa2, sa + sb - hsa2 - 1);


      /* recursively compute a_lo*b_lo into low part of c */
      /* and subtract from T3 */

      KarMul(c, a, hsa, b, hsa, stk);
      KarAdd(T3, c, hsa2 - 1);

      clear(c[hsa2 - 1]);

      /* finally, add T3 * X^{hsa} to c */

      KarAdd(c+hsa, T3, hsa2-1);
   }
   else {
      /* degenerate case */

      GF2X *T;

      T = stk; stk += hsa + sb - 1;

      /* recursively compute b*a_hi into high part of c */

      KarMul(c + hsa, a + hsa, sa - hsa, b, sb, stk);

      /* recursively compute b*a_lo into T */

      KarMul(T, a, hsa, b, sb, stk);

      KarFix(c, T, hsa + sb - 1, hsa);
   }
}

void ExtractBits(_ntl_ulong *cp, const _ntl_ulong *ap, long k, long n)

// extract k bits from a at position n

{
   long sc = (k + NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;

   long wn = n/NTL_BITS_PER_LONG;
   long bn = n - wn*NTL_BITS_PER_LONG;

   long i;

   if (bn == 0) {
      for (i = 0; i < sc; i++)
         cp[i] = ap[i+wn];
   }
   else {
      for (i = 0; i < sc-1; i++)
         cp[i] = (ap[i+wn] >> bn) | (ap[i+wn+1] << (NTL_BITS_PER_LONG - bn));

      if ((k + n) % NTL_BITS_PER_LONG != 0)
         cp[sc-1] = (ap[sc+wn-1] >> bn)|(ap[sc+wn] << (NTL_BITS_PER_LONG - bn));
      else
         cp[sc-1] = ap[sc+wn-1] >> bn;
   }

   long p = k % NTL_BITS_PER_LONG;
   if (p != 0) 
      cp[sc-1] &= ((1UL << p) - 1UL);

}


void KronSubst(GF2X& aa, const GF2EX& a)
{
   long sa = a.rep.length();
   long blocksz = 2*GF2E::degree() - 1;

   long saa = sa*blocksz;

   long wsaa = (saa + NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;

   aa.xrep.SetLength(wsaa+1);

   _ntl_ulong *paa = aa.xrep.elts();


   long i;
   for (i = 0; i < wsaa+1; i++)
      paa[i] = 0;

   for (i = 0; i < sa; i++) 
      ShiftAdd(paa, rep(a.rep[i]).xrep.elts(), rep(a.rep[i]).xrep.length(),
               blocksz*i);

   aa.normalize(); 
}

void KronMul(GF2EX& x, const GF2EX& a, const GF2EX& b)
{
   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());