gf2x.cpp

数值算法库for Windows

#include <NTL/GF2X.h>
#include <NTL/vec_long.h>

#include <NTL/new.h>

NTL_START_IMPL

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


static
long 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;
}

static
istream & 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;   
}    



static
char 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;
}

static
ostream & 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.

inline
void 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;
}

inline
void 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.

static
void 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];
}

static
void 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];