#include <NTL/GF2X.h>#include <NTL/vec_long.h>#include <NTL/new.h>NTL_START_IMPLlong GF2X::HexOutput = 0;void GF2X::SetMaxLength(long n){   if (n < 0) Error("GF2X::SetMaxLength: negative length");   if (n >= (1L << (NTL_BITS_PER_LONG-4)))      Error("GF2X::SetMaxLength: excessive length");   long w = (n + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;   xrep.SetMaxLength(w);}GF2X::GF2X(INIT_SIZE_TYPE, long n){   SetMaxLength(n);}const GF2X& GF2X::zero(){   static GF2X z;   return z;}void GF2X::normalize(){   long n;   const _ntl_ulong *p;   n = xrep.length();   if (n == 0) return;   p = xrep.elts() + (n-1);   while (n > 0 && (*p) == 0) {      p--;      n--;   }   xrep.QuickSetLength(n);}long IsZero(const GF2X& a)    { return a.xrep.length() == 0; }long IsOne(const GF2X& a)   { return a.xrep.length() == 1 && a.xrep[0] == 1; }long IsX(const GF2X& a){   return a.xrep.length() == 1 && a.xrep[0] == 2;}GF2 coeff(const GF2X& a, long i){   if (i < 0) return to_GF2(0);   long wi = i/NTL_BITS_PER_LONG;   if (wi >= a.xrep.length()) return to_GF2(0);   long bi = i - wi*NTL_BITS_PER_LONG;   return to_GF2((a.xrep[wi] & (1UL << bi)) != 0);}GF2 LeadCoeff(const GF2X& a){   if (IsZero(a))      return to_GF2(0);   else      return to_GF2(1);}GF2 ConstTerm(const GF2X& a){   if (IsZero(a))      return to_GF2(0);   else      return to_GF2((a.xrep[0] & 1) != 0);}void set(GF2X& x){   x.xrep.SetLength(1);   x.xrep[0] = 1;}void SetX(GF2X& x){   x.xrep.SetLength(1);   x.xrep[0] = 2;}void SetCoeff(GF2X& x, long i){   if (i < 0) {      Error("SetCoeff: negative index");      return;   }   long n, j;   n = x.xrep.length();   long wi = i/NTL_BITS_PER_LONG;   if (wi >= n) {      x.xrep.SetLength(wi+1);      for (j = n; j <= wi; j++)         x.xrep[j] = 0;   }   long bi = i - wi*NTL_BITS_PER_LONG;   x.xrep[wi] |= (1UL << bi);}   void SetCoeff(GF2X& x, long i, long val){   if (i < 0) {      Error("SetCoeff: negative index");      return;   }   val = val & 1;   if (val) {      SetCoeff(x, i);      return;   }   // we want to clear position i   long n;   n = x.xrep.length();   long wi = i/NTL_BITS_PER_LONG;   if (wi >= n)       return;   long bi = i - wi*NTL_BITS_PER_LONG;   x.xrep[wi] &= ~(1UL << bi);   if (wi == n-1) x.normalize();}void SetCoeff(GF2X& x, long i, GF2 a){   SetCoeff(x, i, rep(a));}void swap(GF2X& a, GF2X& b){   swap(a.xrep, b.xrep);}long deg(const GF2X& aa){   long n = aa.xrep.length();   if (n == 0)      return -1;   _ntl_ulong a = aa.xrep[n-1];   long i = 0;   if (a == 0) Error("GF2X: unnormalized polynomial detected in deg");   while (a>=256)      i += 8, a >>= 8;   if (a >=16)      i += 4, a >>= 4;   if (a >= 4)      i += 2, a >>= 2;   if (a >= 2)      i += 2;   else if (a >= 1)      i++;   return NTL_BITS_PER_LONG*(n-1) + i - 1;}   long operator==(const GF2X& a, const GF2X& b){   return a.xrep == b.xrep;}long operator==(const GF2X& a, long b){   if (b & 1)       return IsOne(a);   else      return IsZero(a);}long operator==(const GF2X& a, GF2 b){   if (b == 1)       return IsOne(a);   else      return IsZero(a);}staticlong FromHex(long c){   if (c >= '0' && c <= '9')      return c - '0';   if (c >= 'A' && c <= 'F')      return 10 + c - 'A';   if (c >= 'a' && c <= 'f')      return 10 + c - 'a';   Error("FromHex: bad arg");   return 0;}staticistream & HexInput(istream& s, GF2X& a){   long n;   long c;   long i;   long val;   GF2X ibuf;   n = 0;   clear(ibuf);   c = s.peek();   while ((c >= '0' && c <= '9') || (c >= 'A' && c <= 'F') ||           (c >= 'a' && c <= 'f'))    {      val = FromHex(c);      for (i = 0; i < 4; i++)         if (val & (1L << i))            SetCoeff(ibuf, n+i);      n += 4;      s.get();      c = s.peek();   }   a = ibuf;   return s;}               istream & operator>>(istream& s, GF2X& a)   {      static ZZ ival;   long c;      if (!s) Error("bad GF2X input");       c = s.peek();     while (c == ' ' || c == '\n' || c == '\t') {        s.get();        c = s.peek();     }     if (c == '0') {      s.get();      c = s.peek();      if (c == 'x' || c == 'X') {         s.get();         return HexInput(s, a);      }      else {         Error("bad GF2X input");      }   }   if (c != '[') {        Error("bad GF2X input");     }     GF2X ibuf;     long n;         n = 0;      clear(ibuf);         s.get();     c = s.peek();     while (c == ' ' || c == '\n' || c == '\t') {        s.get();        c = s.peek();     }     while (c != ']' && c != EOF) {         if (!(s >> ival)) Error("bad GF2X input");      SetCoeff(ibuf, n, to_GF2(ival));      n++;      c = s.peek();        while (c == ' ' || c == '\n' || c == '\t') {           s.get();           c = s.peek();        }     }      if (c == EOF) Error("bad GF2X input");     s.get();       a = ibuf;    return s;   }    staticchar ToHex(long val){   if (val >= 0 && val <= 9)      return char('0' + val);   if (val >= 10 && val <= 15)      return char('a' + val - 10);   Error("ToHex: bad arg");   return 0;}staticostream & HexOutput(ostream& s, const GF2X& a){   s << "0x";   long da = deg(a);   if (da < 0) {      s << '0';      return s;   }   long i, n, val;   val = 0;   n = 0;   for (i = 0; i <= da; i++) {      val = val | (rep(coeff(a, i)) << n);      n++;      if (n == 4) {         s << ToHex(val);         val = 0;         n = 0;      }   }   if (val)       s << ToHex(val);   return s;}ostream& operator<<(ostream& s, const GF2X& a)   {      if (GF2X::HexOutput)      return HexOutput(s, a);   long i, da;      GF2 c;     da = deg(a);      s << '[';         for (i = 0; i <= da; i++) {         c = coeff(a, i);      if (c == 0)         s << "0";      else         s << "1";      if (i < da) s << " ";      }         s << ']';            return s;   }   void random(GF2X& x, long n){   if (n < 0) Error("GF2X random: negative length");   if (n >= (1L << (NTL_BITS_PER_LONG-4)))      Error("GF2X random: excessive length");   long wl = (n+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;   x.xrep.SetLength(wl);   long i;   for (i = 0; i < wl-1; i++) {      x.xrep[i] = RandomWord();   }   if (n > 0) {      long pos = n % NTL_BITS_PER_LONG;      if (pos == 0) pos = NTL_BITS_PER_LONG;      x.xrep[wl-1] = RandomBits_ulong(pos);   }   x.normalize();}void add(GF2X& x, const GF2X& a, const GF2X& b){   long sa = a.xrep.length();   long sb = b.xrep.length();   long i;   if (sa == sb) {      x.xrep.SetLength(sa);      if (sa == 0) return;      _ntl_ulong *xp = x.xrep.elts();      const _ntl_ulong *ap = a.xrep.elts();      const _ntl_ulong *bp = b.xrep.elts();      for (i = 0; i < sa; i++)         xp[i] = ap[i] ^ bp[i];      i = sa-1;      while (i >= 0 && !xp[i]) i--;      x.xrep.QuickSetLength(i+1);   }      else if (sa < sb) {      x.xrep.SetLength(sb);      _ntl_ulong *xp = x.xrep.elts();      const _ntl_ulong *ap = a.xrep.elts();      const _ntl_ulong *bp = b.xrep.elts();      for (i = 0; i < sa; i++)         xp[i] = ap[i] ^ bp[i];      for (; i < sb; i++)         xp[i] = bp[i];   }   else { // sa > sb      x.xrep.SetLength(sa);      _ntl_ulong *xp = x.xrep.elts();      const _ntl_ulong *ap = a.xrep.elts();      const _ntl_ulong *bp = b.xrep.elts();      for (i = 0; i < sb; i++)         xp[i] = ap[i] ^ bp[i];      for (; i < sa; i++)         xp[i] = ap[i];   }}// This mul1 routine (for 32 x 32 bit multiplies) // was arrived at after a good deal of experimentation.// Several people have advocated that the "best" way// is to reduce it via karastuba to 9 8x8 multiplies, implementing// the latter via table look-up (a huge table).  Although such a mul1// would indeed have fewer machine instructions, my// experiments on a PowerPC and on a Pentium indicate// that the memory accesses slow it down.  On these machines// the following mul1 is faster by a factor of about 1.5.// inlining mul1 seems to help with g++; morever,// the IBM xlC compiler inlines it anyway.inlinevoid mul1(_ntl_ulong *c, _ntl_ulong a, _ntl_ulong b){   _ntl_ulong hi, lo;   _ntl_ulong A[4];   A[0] = 0;   A[1] = a & ((1UL << (NTL_BITS_PER_LONG-1))-1UL);   A[2] = a << 1;   A[3] = A[1] ^ A[2];   lo = A[b>>(NTL_BITS_PER_LONG-2)];   hi = lo >> (NTL_BITS_PER_LONG-2);    lo = (lo << 2) ^ A[(b >> (NTL_BITS_PER_LONG-4)) & 3];   // The following code is included from mach_desc.h   // and handles *any* word size   NTL_BB_MUL_CODE   hi = (hi << 2)|(lo >> (NTL_BITS_PER_LONG-2));    lo = (lo << 2) ^ A[b & 3];   if (a >> (NTL_BITS_PER_LONG-1)) {      hi = hi ^ (b >> 1);      lo = lo ^ (b << (NTL_BITS_PER_LONG-1));   }   c[0] = lo;   c[1] = hi;}inlinevoid mul_half(_ntl_ulong *c, _ntl_ulong a, _ntl_ulong b){   _ntl_ulong hi, lo;   _ntl_ulong A[4];   A[0] = 0;   A[1] = a & ((1UL << (NTL_BITS_PER_LONG-1))-1UL);   A[2] = a << 1;   A[3] = A[1] ^ A[2];   lo = A[b>>(NTL_BITS_PER_LONG/2-2)];   hi = lo >> (NTL_BITS_PER_LONG-2);    lo = (lo << 2) ^ A[(b >> (NTL_BITS_PER_LONG/2-4)) & 3];   NTL_BB_HALF_MUL_CODE   hi = (hi << 2)|(lo >> (NTL_BITS_PER_LONG-2));    lo = (lo << 2) ^ A[b & 3];   if (a >> (NTL_BITS_PER_LONG-1)) {      hi = hi ^ (b >> 1);      lo = lo ^ (b << (NTL_BITS_PER_LONG-1));   }   c[0] = lo;   c[1] = hi;}// mul2...mul8 hard-code 2x2...8x8 word multiplies.// I adapted these routines from LiDIA.// Inlining these seems to hurt, not help.staticvoid mul2(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b){   _ntl_ulong hs0, hs1;   _ntl_ulong hl2[2];   hs0 = a[0] ^ a[1];   hs1 = b[0] ^ b[1];   mul1(c, a[0], b[0]);   mul1(c+2, a[1], b[1]);   mul1(hl2, hs0, hs1);   hl2[0] = hl2[0] ^ c[0] ^ c[2];   hl2[1] = hl2[1] ^ c[1] ^ c[3];   c[1] ^= hl2[0];   c[2] ^= hl2[1];}staticvoid mul3 (_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b){   _ntl_ulong hs0[2], hs1[2];   _ntl_ulong hl2[4];   hs0[0] = a[0] ^ a[2];    hs0[1] = a[1];   hs1[0] = b[0] ^ b[2];    hs1[1] = b[1];   mul2(c, a, b);    mul1(c+4, a[2], b[2]);   mul2(hl2, hs0, hs1);    hl2[0] = hl2[0] ^ c[0] ^ c[4];    hl2[1] = hl2[1] ^ c[1] ^ c[5];