gf2x1.cpp

一个比较通用的大数运算库

const long GF2X_DIV_CROSS = 100; 

void DivRem(GF2X& q, GF2X& r, const GF2X& a, const GF2X& b)
{
   long sa = a.xrep.length();
   long sb = b.xrep.length();

   if (sb < GF2X_DIV_CROSS || sa-sb < GF2X_DIV_CROSS)
      PlainDivRem(q, r, a, b);
   else if (sa < 4*sb)
      UseMulDivRem(q, r, a, b);
   else {
      GF2XModulus B;
      build(B, b);
      DivRem(q, r, a, B);
   }
}

void div(GF2X& q, const GF2X& a, const GF2X& b)
{
   long sa = a.xrep.length();
   long sb = b.xrep.length();

   if (sb < GF2X_DIV_CROSS || sa-sb < GF2X_DIV_CROSS)
      PlainDiv(q, a, b);
   else if (sa < 4*sb)
      UseMulDiv(q, a, b);
   else {
      GF2XModulus B;
      build(B, b);
      div(q, a, B);
   }
}

void rem(GF2X& r, const GF2X& a, const GF2X& b)
{
   long sa = a.xrep.length();
   long sb = b.xrep.length();

   if (sb < GF2X_DIV_CROSS || sa-sb < GF2X_DIV_CROSS)
      PlainRem(r, a, b);
   else if (sa < 4*sb)
      UseMulRem(r, a, b);
   else {
      GF2XModulus B;
      build(B, b);
      rem(r, a, B);
   }
}


static inline 
void swap(_ntl_ulong_ptr& a, _ntl_ulong_ptr& b)  
{  _ntl_ulong_ptr t;  t = a; a = b; b = t; }




static
void BaseGCD(GF2X& d, const GF2X& a_in, const GF2X& b_in)
{
   GF2XRegister(a);
   GF2XRegister(b);
   
   if (IsZero(a_in)) {
      d = b_in;
      return;
   }

   if (IsZero(b_in)) {
      d = a_in;
      return;
   }
      
   a.xrep.SetMaxLength(a_in.xrep.length()+1);
   b.xrep.SetMaxLength(b_in.xrep.length()+1);

   a = a_in;
   b = b_in;

   _ntl_ulong *ap = a.xrep.elts();
   _ntl_ulong *bp = b.xrep.elts();

   long da = deg(a);
   long wa = da/NTL_BITS_PER_LONG;
   long ba = da - wa*NTL_BITS_PER_LONG;

   long db = deg(b);
   long wb = db/NTL_BITS_PER_LONG;
   long bb = db - wb*NTL_BITS_PER_LONG;

   long parity = 0;

   for (;;) {
      if (da < db) {
         swap(ap, bp);
         swap(da, db);
         swap(wa, wb);
         swap(ba, bb);
         parity = 1 - parity;
      }

      // da >= db

      if (db == -1) break;

      ShiftAdd(ap, bp, wb+1, da-db);

      _ntl_ulong msk = 1UL << ba;
      _ntl_ulong aa = ap[wa];

      while ((aa & msk) == 0) {
         da--;
         msk = msk >> 1;
         ba--;
         if (!msk) {
            wa--;
            ba = NTL_BITS_PER_LONG-1;
            msk = 1UL << (NTL_BITS_PER_LONG-1);
            if (wa < 0) break;
            aa = ap[wa];
         }
      }
   }

   a.normalize();
   b.normalize();

   if (!parity) {
      d = a;
   }
   else {
      d = b;
   }
}


void GCD(GF2X& d, const GF2X& a, const GF2X& b)
{
   long sa = a.xrep.length();
   long sb = b.xrep.length();

   if (sb >= 10 && 2*sa > 3*sb) {
      GF2XRegister(r);

      rem(r, a, b);
      BaseGCD(d, b, r);
   }
   else if (sa >= 10 && 2*sb > 3*sa) {
      GF2XRegister(r);

      rem(r, b, a);
      BaseGCD(d, a, r);
   }
   else {
      BaseGCD(d, a, b);
   }
}





#define XX_STEP(ap,da,wa,ba,rp,sr,bp,db,wb,bb,sp,ss)  \
      long delta = da-db;  \
  \
      if (delta == 0) {  \
         long i;  \
         for (i = wb; i >= 0; i--) ap[i] ^= bp[i];  \
         for (i = ss-1; i >= 0; i--) rp[i] ^= sp[i];  \
         if (ss > sr) sr = ss; \
      }  \
      else if (delta == 1) {  \
         long i; \
         _ntl_ulong tt, tt1;  \
  \
         tt = bp[wb] >> (NTL_BITS_PER_LONG-1);  \
         if (tt) ap[wb+1] ^= tt;  \
         tt = bp[wb];  \
         for (i = wb; i >= 1; i--)  \
            tt1 = bp[i-1], ap[i] ^= (tt << 1) | (tt1 >> (NTL_BITS_PER_LONG-1)),  \
            tt = tt1; \
         ap[0] ^= tt << 1;  \
  \
         if (ss > 0) {  \
            long t = ss; \
            tt = sp[ss-1] >> (NTL_BITS_PER_LONG-1);  \
            if (tt) rp[ss] ^= tt, t++;  \
            tt = sp[ss-1]; \
            for (i = ss-1; i >= 1; i--)  \
               tt1=sp[i-1],  \
               rp[i] ^= (tt << 1) | (tt1 >> (NTL_BITS_PER_LONG-1)),  \
               tt = tt1; \
            rp[0] ^= tt << 1;  \
            if (t > sr) sr = t; \
         }  \
      }  \
      else if (delta < NTL_BITS_PER_LONG) {  \
         long i; \
         _ntl_ulong tt, tt1;  \
         long rdelta = NTL_BITS_PER_LONG-delta; \
  \
         tt = bp[wb] >> rdelta;  \
         if (tt) ap[wb+1] ^= tt;  \
         tt=bp[wb]; \
         for (i = wb; i >= 1; i--)  \
            tt1=bp[i-1], ap[i] ^= (tt << delta) | (tt1 >> rdelta),  \
            tt=tt1; \
         ap[0] ^= tt << delta;  \
  \
         if (ss > 0) {  \
            long t = ss; \
            tt = sp[ss-1] >> rdelta;  \
            if (tt) rp[ss] ^= tt, t++;  \
            tt=sp[ss-1]; \
            for (i = ss-1; i >= 1; i--)  \
               tt1=sp[i-1], rp[i] ^= (tt << delta) | (tt1 >> rdelta),  \
               tt=tt1; \
            rp[0] ^= tt << delta;  \
            if (t > sr) sr = t; \
         }  \
      }  \
      else {  \
         ShiftAdd(ap, bp, wb+1, da-db);  \
         ShiftAdd(rp, sp, ss, da-db);  \
         long t = ss + (da-db+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;  \
         if (t > sr) {  \
            while (t > 0 && rp[t-1] == 0) t--;   \
            sr = t;  \
         }  \
      } \
  \
      _ntl_ulong msk = 1UL << ba;  \
      _ntl_ulong aa = ap[wa];  \
  \
      while ((aa & msk) == 0) {  \
         da--;  \
         msk = msk >> 1;  \
         ba--;  \
         if (!msk) {  \
            wa--;  \
            ba = NTL_BITS_PER_LONG-1;  \
            msk = 1UL << (NTL_BITS_PER_LONG-1);  \
            if (wa < 0) break;  \
            aa = ap[wa];  \
         }  \
      }  \




static
void XXGCD(GF2X& d, GF2X& r_out, const GF2X& a_in, const GF2X& b_in)
{
   GF2XRegister(a);
   GF2XRegister(b);
   GF2XRegister(r);
   GF2XRegister(s);

   if (IsZero(b_in)) {
      d = a_in;
      set(r_out);
      return;
   }

   if (IsZero(a_in)) {
      d = b_in;
      clear(r_out);
      return;
   }
      
   a.xrep.SetMaxLength(a_in.xrep.length()+1);
   b.xrep.SetMaxLength(b_in.xrep.length()+1);

   long max_sz = max(a_in.xrep.length(), b_in.xrep.length());
   r.xrep.SetLength(max_sz+1);
   s.xrep.SetLength(max_sz+1);

   _ntl_ulong *rp = r.xrep.elts();
   _ntl_ulong *sp = s.xrep.elts();

   long i;
   for (i = 0; i <= max_sz; i++) {
      rp[i] = sp[i] = 0;
   }

   rp[0] = 1;

   long sr = 1;
   long ss = 0;

   a = a_in;
   b = b_in;

   _ntl_ulong *ap = a.xrep.elts();
   _ntl_ulong *bp = b.xrep.elts();

   long da = deg(a);
   long wa = da/NTL_BITS_PER_LONG;
   long ba = da - wa*NTL_BITS_PER_LONG;

   long db = deg(b);
   long wb = db/NTL_BITS_PER_LONG;
   long bb = db - wb*NTL_BITS_PER_LONG;

   long parity = 0;


   for (;;) {
      if (da == -1 || db == -1) break;

      if (da < db || (da == db && parity)) {
         if (da < db && !parity) parity = 1;
         XX_STEP(bp,db,wb,bb,sp,ss,ap,da,wa,ba,rp,sr)

      }
      else {
         parity = 0;
         XX_STEP(ap,da,wa,ba,rp,sr,bp,db,wb,bb,sp,ss)
      }
   }

   a.normalize();
   b.normalize();
   r.normalize();
   s.normalize();

   if (db == -1) {
      d = a;
      r_out = r;
   }
   else {
      d = b;
      r_out = s;
   }
}



static
void BaseXGCD(GF2X& d, GF2X& s, GF2X& t, const GF2X& a, const GF2X& b)
{
   if (IsZero(b)) {
      d = a;
      set(s);
      clear(t);
   }
   else {
      GF2XRegister(t1);
      GF2XRegister(b1);

      b1 = b;
      XXGCD(d, s, a, b);
      mul(t1, a, s);
      add(t1, t1, d);
      div(t, t1, b1);
   }
}




void XGCD(GF2X& d, GF2X& s, GF2X& t, const GF2X& a, const GF2X& b)
{
   long sa = a.xrep.length();
   long sb = b.xrep.length();


   if (sb >= 10 && 2*sa > 3*sb) {
      GF2XRegister(r);
      GF2XRegister(q);
      GF2XRegister(s1);
      GF2XRegister(t1);


      DivRem(q, r, a, b);
      BaseXGCD(d, s1, t1, b, r);

      
      mul(r, t1, q);
      add(r, r, s1);  // r = s1 - t1*q, but sign doesn't matter

      s = t1;
      t = r;   
   }
   else if (sa >= 10 && 2*sb > 3*sa) {
      GF2XRegister(r);
      GF2XRegister(q);
      GF2XRegister(s1);
      GF2XRegister(t1);


      DivRem(q, r, b, a);
      BaseXGCD(d, s1, t1, a, r);

      
      mul(r, t1, q);
      add(r, r, s1);  // r = s1 - t1*q, but sign doesn't matter

      t = t1;
      s = r;  
   }
   else {
      BaseXGCD(d, s, t, a, b);
   }

}




static
void BaseInvMod(GF2X& d, GF2X& s, const GF2X& a, const GF2X& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args");

   long sa = a.xrep.length();
   long sf = f.xrep.length();

   if (sa >= 10 && 2*sf > 3*sa) {
      GF2XRegister(t);

      XGCD(d, s, t, a, f);
   }
   else {
      XXGCD(d, s, a, f);
   }

}



void InvMod(GF2X& c, const GF2X& a, const GF2X& f)
{ 
   GF2XRegister(d);
   GF2XRegister(s);
   BaseInvMod(d, s, a, f);

   if (!IsOne(d)) Error("InvMod: inverse undefined");

   c = s;
}



long InvModStatus(GF2X& c, const GF2X& a, const GF2X& f)
{ 
   GF2XRegister(d);
   GF2XRegister(s);
   BaseInvMod(d, s, a, f);

   if (!IsOne(d)) {
      c = d;
      return 1;
   }

   c = s;
   return 0;
}



   
void diff(GF2X& c, const GF2X& a)
{
   RightShift(c, a, 1);
   
   // clear odd coeffs

   long dc = deg(c);
   long i;
   for (i = 1; i <= dc; i += 2)
      SetCoeff(c, i, 0);
}

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

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

void conv(GF2X& x, const vec_GF2& a)
{
   x.xrep = a.rep;
   x.normalize();
}

void conv(vec_GF2& x, const GF2X& a)
{
   VectorCopy(x, a, deg(a)+1);
}

void VectorCopy(vec_GF2& x, const GF2X& a, long n)
{
   if (n < 0) Error("VectorCopy: negative length"); 

   if (NTL_OVERFLOW(n, 1, 0))
      Error("overflow in VectorCopy");

   long wa = a.xrep.length();
   long wx = (n + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;

   long wmin = min(wa, wx);

   x.SetLength(n);

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

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

   if (wa < wx) {
      for (i = wa; i < wx; i++)
         xp[i] = 0;
   }
   else {
      long p = n % NTL_BITS_PER_LONG;
      if (p != 0)
         xp[wx-1] &= (1UL << p) - 1UL;
   }
}


void add(GF2X& c, const GF2X& a, long b)
{
   c = a;
   if (b & 1) {
      long n = c.xrep.length();
      if (n == 0) 
         set(c);
      else {
         c.xrep[0] ^= 1;
         if (n == 1 && !c.xrep[0]) c.xrep.SetLength(0);
      }
   }
}

void add(GF2X& c, const GF2X& a, GF2 b)
{
   add(c, a, rep(b));
}


void MulTrunc(GF2X& c, const GF2X& a, const GF2X& b, long n)
{
   GF2XRegister(t);

   mul(t, a, b);
   trunc(c, t, n);
}

void SqrTrunc(GF2X& c, const GF2X& a, long n)
{
   GF2XRegister(t);

   sqr(t, a);
   trunc(c, t, n);
}


long divide(GF2X& q, const GF2X& a, const GF2X& b)
{
   if (IsZero(b)) {
      if (IsZero(a)) {
         clear(q);
         return 1;
      }
      else
         return 0;
   }

   GF2XRegister(lq);
   GF2XRegister(r);

   DivRem(lq, r, a, b);
   if (!IsZero(r)) return 0;
   q = lq;
   return 1;
}

long divide(const GF2X& a, const GF2X& b)
{
   if (IsZero(b)) return IsZero(a);
   GF2XRegister(r);
   rem(r, a, b);
   if (!IsZero(r)) return 0;
   return 1;
}



/*** modular composition routines and data structures ***/


void InnerProduct(GF2X& x, const GF2X& v, long dv, long low, long high, 
                   const vec_GF2X& H, long n, WordVector& t)
{
   long i, j;

   _ntl_ulong *tp = t.elts();

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


   long w_low = low/NTL_BITS_PER_LONG;
   long b_low = low - w_low*NTL_BITS_PER_LONG;

   
   const _ntl_ulong *vp = &v.xrep[w_low];
   _ntl_ulong msk = 1UL << b_low;
   _ntl_ulong vv = *vp;

   high = min(high, dv);

   i = low;
   for (;;) {
      if (vv & msk) {
         const WordVector& h = H[i-low].xrep;
         long m = h.length();
         const _ntl_ulong *hp = h.elts();
         for (j = 0; j < m; j++)
            tp[j] ^= hp[j];
      }

      i++;
      if (i > high) break;

      msk = msk << 1;
      if (!msk) {
         msk = 1UL;