gf2x1.c

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

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

#include <NTL/new.h>

NTL_START_IMPL

/********** data structures for accesss to GF2XRegisters ************/

static GF2X GF2XRegisterVec[32];
static long GF2XRegisterTop = 0;


class GF2XRegisterType {
public:

GF2X *xrep;

GF2XRegisterType()
{ xrep = &GF2XRegisterVec[GF2XRegisterTop]; GF2XRegisterTop++; }

~GF2XRegisterType()
{ GF2XRegisterTop--; }

operator GF2X& () { return *xrep; }

};

#define GF2XRegister(a) GF2XRegisterType GF2XReg__ ## a ; GF2X& a = GF2XReg__ ## a







static vec_GF2X stab;  // used by PlainDivRem and PlainRem

static WordVector GF2X_rembuf;


void PlainDivRem(GF2X& q, GF2X& r, const GF2X& a, const GF2X& b)
{
   long da, sa, posa, db, sb, posb, dq, sq, posq;

   da = deg(a);
   db = deg(b);

   if (db < 0) Error("GF2X: division by zero");

   if (da < db) {
      r = a;
      clear(q);
      return;
   }

   sa = a.xrep.length();
   posa = da - NTL_BITS_PER_LONG*(sa-1);
   sb = b.xrep.length();
   posb = db - NTL_BITS_PER_LONG*(sb-1);

   dq = da - db;
   sq = dq/NTL_BITS_PER_LONG + 1;
   posq = dq - NTL_BITS_PER_LONG*(sq-1);

   _ntl_ulong *ap;
   if (&r == &a)
      ap = r.xrep.elts();
   else {
      GF2X_rembuf = a.xrep;
      ap = GF2X_rembuf.elts();
   }

   stab.SetLength(NTL_BITS_PER_LONG);
   long i;

   stab[posb] = b;
   for (i = 1; i <= min(dq, NTL_BITS_PER_LONG-1); i++) 
      MulByX(stab[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG], 
             stab[((_ntl_ulong)(posb+i-1))%NTL_BITS_PER_LONG]);

   _ntl_ulong *stab_ptr[NTL_BITS_PER_LONG];
   long stab_cnt[NTL_BITS_PER_LONG];

   for (i = 0; i <= min(dq, NTL_BITS_PER_LONG-1); i++) {
      WordVector& st = stab[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG].xrep;
      long k = st.length();
      stab_ptr[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG] = &st[k-1];
      stab_cnt[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG] = -k+1;
   }

   q.xrep.SetLength(sq);
   _ntl_ulong *qp = q.xrep.elts();
   for (i = 0; i < sq; i++)
      qp[i] = 0;

   _ntl_ulong *atop = &ap[sa-1];
   _ntl_ulong *qtop = &qp[sq-1];
   _ntl_ulong *stab_top;

   while (da >= db) {
      if (atop[0] & (1UL << posa)) {
         qtop[0] |= (1UL << posq);
         stab_top = stab_ptr[posa];
         for (i = stab_cnt[posa]; i <= 0; i++)
            atop[i] ^= stab_top[i];
      }

      da--;
      posa--;
      if (posa < 0) {
         posa = NTL_BITS_PER_LONG-1;
         atop--;
      }

      posq--;
      if (posq < 0) {
         posq = NTL_BITS_PER_LONG-1;
         qtop--;
      }
   }

   if (posb == 0) sb--;

   r.xrep.SetLength(sb);
   if (&r != &a) {
      _ntl_ulong *rp = r.xrep.elts();
      for (i = 0; i < sb; i++)
         rp[i] = ap[i];
   }
   r.normalize();
}



void PlainDiv(GF2X& q, const GF2X& a, const GF2X& b)
{
   GF2XRegister(r);
   PlainDivRem(q, r, a, b);
}


void PlainRem(GF2X& r, const GF2X& a, const GF2X& b)
{
   long da, sa, posa, db, sb, posb;

   da = deg(a);
   db = deg(b);

   if (db < 0) Error("GF2X: division by zero");

   if (da < db) {
      r = a;
      return;
   }

   sa = a.xrep.length();
   posa = da - NTL_BITS_PER_LONG*(sa-1);
   sb = b.xrep.length();
   posb = db - NTL_BITS_PER_LONG*(sb-1);

   _ntl_ulong *ap;
   if (&r == &a)
      ap = r.xrep.elts();
   else {
      GF2X_rembuf = a.xrep;
      ap = GF2X_rembuf.elts();
   }

   stab.SetLength(NTL_BITS_PER_LONG);
   long i;

   stab[posb] = b;
   for (i = 1; i <= min(da-db, NTL_BITS_PER_LONG-1); i++) 
      MulByX(stab[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG], 
             stab[((_ntl_ulong)(posb+i-1))%NTL_BITS_PER_LONG]);

   _ntl_ulong *stab_ptr[NTL_BITS_PER_LONG];
   long stab_cnt[NTL_BITS_PER_LONG];

   for (i = 0; i <= min(da-db, NTL_BITS_PER_LONG-1); i++) {
      WordVector& st = stab[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG].xrep;
      long k = st.length();
      stab_ptr[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG] = &st[k-1];
      stab_cnt[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG] = -k+1;
   }


   _ntl_ulong *atop = &ap[sa-1];
   _ntl_ulong *stab_top;

   while (da >= db) {
      if (atop[0] & (1UL << posa)) {
         stab_top = stab_ptr[posa];
         for (i = stab_cnt[posa]; i <= 0; i++)
            atop[i] ^= stab_top[i];
      }

      da--;
      posa--;
      if (posa < 0) {
         posa = NTL_BITS_PER_LONG-1;
         atop--;
      }
   }

   if (posb == 0) sb--;

   r.xrep.SetLength(sb);
   if (&r != &a) {
      _ntl_ulong *rp = r.xrep.elts();
      for (i = 0; i < sb; i++)
         rp[i] = ap[i];
   }
   r.normalize();
}

#define MASK8 ((1UL << 8)-1UL)

static _ntl_ulong invtab[128] = {
1UL, 255UL, 85UL, 219UL, 73UL, 151UL, 157UL, 51UL, 17UL, 175UL,
69UL, 139UL, 89UL, 199UL, 141UL, 99UL, 33UL, 95UL, 117UL, 123UL,
105UL, 55UL, 189UL, 147UL, 49UL, 15UL, 101UL, 43UL, 121UL, 103UL,
173UL, 195UL, 65UL, 191UL, 21UL, 155UL, 9UL, 215UL, 221UL, 115UL,
81UL, 239UL, 5UL, 203UL, 25UL, 135UL, 205UL, 35UL, 97UL, 31UL,
53UL, 59UL, 41UL, 119UL, 253UL, 211UL, 113UL, 79UL, 37UL, 107UL,
57UL, 39UL, 237UL, 131UL, 129UL, 127UL, 213UL, 91UL, 201UL, 23UL,
29UL, 179UL, 145UL, 47UL, 197UL, 11UL, 217UL, 71UL, 13UL, 227UL,
161UL, 223UL, 245UL, 251UL, 233UL, 183UL, 61UL, 19UL, 177UL, 143UL,
229UL, 171UL, 249UL, 231UL, 45UL, 67UL, 193UL, 63UL, 149UL, 27UL,
137UL, 87UL, 93UL, 243UL, 209UL, 111UL, 133UL, 75UL, 153UL, 7UL,
77UL, 163UL, 225UL, 159UL, 181UL, 187UL, 169UL, 247UL, 125UL, 83UL,
241UL, 207UL, 165UL, 235UL, 185UL, 167UL, 109UL, 3UL };



void NewtonInvTrunc(GF2X& c, const GF2X& a, long e)
{
   if (e == 1) {
      set(c);
      return;
   }

   static vec_long E;
   E.SetLength(0);
   append(E, e);
   while (e > 8) {
      e = (e+1)/2;
      append(E, e);
   }

   long L = E.length();

   GF2XRegister(g);
   GF2XRegister(g0);
   GF2XRegister(g1);
   GF2XRegister(g2);

   g.xrep.SetMaxLength((e+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG + 1);
   g0.xrep.SetMaxLength((e+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG + 1);
   g1.xrep.SetMaxLength(((3*e+1)/2+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG + 1);
   g2.xrep.SetMaxLength((e+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG + 1);

   g.xrep.SetLength(1);
   g.xrep[0] = invtab[(a.xrep[0] & MASK8) >> 1] & ((1UL<<e)-1UL);

   long i;

   for (i = L-1; i > 0; i--) {
      // lift from E[i] to E[i-1]

      long k = E[i];
      long l = E[i-1]-E[i];

      trunc(g0, a, k+l);

      mul(g1, g0, g);
      RightShift(g1, g1, k);
      trunc(g1, g1, l);

      mul(g2, g1, g);
      trunc(g2, g2, l);
      LeftShift(g2, g2, k);

      add(g, g, g2);
   }

   c = g;
}

void InvTrunc(GF2X& c, const GF2X& a, long e)
{
   if (ConstTerm(a) == 0 || e < 0)
      Error("inv: bad args");

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

   if (e == 0) {
      clear(c);
      return;
   }

   NewtonInvTrunc(c, a, e);
}



static long weight1(_ntl_ulong a)
{
   long res = 0;
   while (a) {
      if (a & 1) res ++;
      a >>= 1;
   }
   return res;
}

long weight(const GF2X& a)
{
   long wlen = a.xrep.length();
   long res = 0;
   long i;
   for (i = 0; i < wlen; i++)
      res += weight1(a.xrep[i]);

   return res;
}



static
void SparsityCheck(const GF2X& f, long& k3, long& k2, long& k1)
{
   long w = weight(f);
   if (w != 3 && w != 5) {
      k3 = 0;
      return;
   }

   if (ConstTerm(f) != 1) {
      k3 = 0;
      return;
   }

   GF2X g = f;

   long n = deg(f);

   trunc(g, g, n);
   
   long t = deg(g);

   if (n-t < NTL_BITS_PER_LONG || t > (n+1)/2) {
      k3 = 0;
      return;
   }

   if (w == 3) {
      k3 = t;
      k2 = 0;
      return;
   }

   k3 = t;
   trunc(g, g, t);
   t = deg(g);
   k2 = t;
   trunc(g, g, t);
   t = deg(g);
   k1 = t;
}




const long GF2X_MOD_PLAIN = 0;
const long GF2X_MOD_MUL = 1;
const long GF2X_MOD_SPECIAL = 2;
const long GF2X_MOD_TRI = 3;
const long GF2X_MOD_PENT = 4;

void build(GF2XModulus& F, const GF2X& f)
{
   long n = deg(f);
   long i;

   if (n <= 0) Error("build(GF2XModulus,GF2X): deg(f) <= 0");

   F.tracevec.SetLength(0);

   F.f = f;
   F.n = n;
   F.sn = f.xrep.length();

   long sb = F.sn;
   long posb = n - NTL_BITS_PER_LONG*(sb-1);

   F.posn = posb; 

   if (F.posn > 0) {
      F.size = F.sn;
      F.msk = (1UL << F.posn) - 1UL;
   }
   else {
      F.size = F.sn-1;
      F.msk = ~0UL;
   }

   SparsityCheck(f, F.k3, F.k2, F.k1);

   if (F.k3 != 0) {
      if (F.k2 == 0)
         F.method = GF2X_MOD_TRI;
      else
         F.method = GF2X_MOD_PENT;

      return;
   }


   GF2X f0;
   trunc(f0, f, n);
   long deg_f0 = deg(f0);

   if (F.sn > 1 && deg_f0 < NTL_BITS_PER_LONG 
       && deg_f0 >= NTL_BITS_PER_LONG/2) {
      if (F.size >= 6)
         F.method = GF2X_MOD_MUL;
      else
         F.method = GF2X_MOD_SPECIAL;
   }
   else if (F.sn > 1 && deg_f0 < NTL_BITS_PER_LONG/2) {
      if (F.size >= 4)
         F.method = GF2X_MOD_MUL;
      else
         F.method = GF2X_MOD_SPECIAL;
   }
   else if (F.size >= 8)
      F.method = GF2X_MOD_MUL;
   else 
      F.method = GF2X_MOD_PLAIN;
      

   if (F.method == GF2X_MOD_SPECIAL) {
      if (!F.stab_cnt) F.stab_cnt = NTL_NEW_OP long[NTL_BITS_PER_LONG];
      long *stab_cnt = F.stab_cnt;
      if (!stab_cnt) Error("out of memory");

      if (!F.stab1) F.stab1 = NTL_NEW_OP _ntl_ulong[2*NTL_BITS_PER_LONG];
      _ntl_ulong *stab1 = F.stab1;
      if (!stab1) Error("out of memory");

      stab1[posb<<1] = f.xrep[0];
      stab1[(posb<<1)+1] = 0;

      stab_cnt[posb] = -sb+1;

      for (i = 1; i < NTL_BITS_PER_LONG; i++) {
         long kk0 = ((_ntl_ulong)(posb+i-1))%NTL_BITS_PER_LONG;
         long kk1 = ((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG;

         stab1[kk1<<1] = stab1[kk0<<1] << 1;
         stab1[(kk1<<1)+1] = (stab1[(kk0<<1)+1] << 1) 
                          | (stab1[kk0<<1] >> (NTL_BITS_PER_LONG-1));

         if (kk1 < posb) 
            stab_cnt[kk1] = -sb;
         else
            stab_cnt[kk1] = -sb+1;
      }
   }
   else if (F.method == GF2X_MOD_PLAIN) {
      vec_GF2X& stab = F.stab;
      stab.SetLength(NTL_BITS_PER_LONG);


      if (!F.stab_ptr) F.stab_ptr = NTL_NEW_OP _ntl_ulong_ptr[NTL_BITS_PER_LONG];
      _ntl_ulong **stab_ptr = F.stab_ptr;
      if (!stab_ptr) Error("out of memory");

      if (!F.stab_cnt) F.stab_cnt = NTL_NEW_OP long[NTL_BITS_PER_LONG];
      long *stab_cnt = F.stab_cnt;
      if (!stab_cnt) Error("out of memory");
      
   
      stab[posb] = f;
      for (i = 1; i < NTL_BITS_PER_LONG; i++) 
         MulByX(stab[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG], 
                stab[((_ntl_ulong)(posb+i-1))%NTL_BITS_PER_LONG]);
   
   
      for (i = 0; i < NTL_BITS_PER_LONG; i++) {
         WordVector& st = stab[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG].xrep;
         long k = st.length();
         stab_ptr[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG] = &st[k-1];
         stab_cnt[((_ntl_ulong)(posb+i))%NTL_BITS_PER_LONG] = -k+1;
      }
   }
   else if (F.method == GF2X_MOD_MUL) {
      GF2X P1, P2;

      CopyReverse(P1, f, n);
      InvTrunc(P2, P1, n-1);
      CopyReverse(P1, P2, n-2);
      trunc(F.h0, P1, n-2);
      F.f0 = f0;
   }
}

GF2XModulus::GF2XModulus()
{
   n = -1;
   method = GF2X_MOD_PLAIN;
   stab_ptr = 0;
   stab_cnt = 0;
   stab1 = 0;
}


// The following two routines are total spaghetti...unfortunately,
// cleaning them up would require too much re-coding in other
// places.

GF2XModulus::GF2XModulus(const GF2XModulus& F) :
   f(F.f), n(F.n), sn(F.sn), posn(F.posn), k3(F.k3), k2(F.k2), k1(F.k1),
   size(F.size), 
   msk(F.msk), method(F.method), stab(F.stab), h0(F.h0), f0(F.f0),
   stab_cnt(0), stab_ptr(0), stab1(0), tracevec(F.tracevec)
{
   if (method == GF2X_MOD_SPECIAL) {
      long i;
      stab1 = NTL_NEW_OP _ntl_ulong[2*NTL_BITS_PER_LONG];
      if (!stab1) Error("GF2XModulus: out of memory");
      for (i = 0; i < 2*NTL_BITS_PER_LONG; i++)
         stab1[i] = F.stab1[i];
      stab_cnt = NTL_NEW_OP long[NTL_BITS_PER_LONG];
      if (!stab_cnt) Error("GF2XModulus: out of memory");
      for (i = 0; i < NTL_BITS_PER_LONG; i++)
         stab_cnt[i] = F.stab_cnt[i];
   }
   else if (method == GF2X_MOD_PLAIN) {
      long i;

      if (F.stab_cnt) {
         stab_cnt = NTL_NEW_OP long[NTL_BITS_PER_LONG];
         if (!stab_cnt) Error("GF2XModulus: out of memory");
         for (i = 0; i < NTL_BITS_PER_LONG; i++)
            stab_cnt[i] = F.stab_cnt[i];
      }

      if (F.stab_ptr) {
         stab_ptr = NTL_NEW_OP _ntl_ulong_ptr[NTL_BITS_PER_LONG];
         if (!stab_ptr) Error("GF2XModulus: out of memory");
      
         for (i = 0; i < NTL_BITS_PER_LONG; i++) {
            WordVector& st = stab[((_ntl_ulong)(posn+i))%NTL_BITS_PER_LONG].xrep;
            long k = st.length();
            stab_ptr[((_ntl_ulong)(posn+i))%NTL_BITS_PER_LONG] = &st[k-1];
            stab_cnt[((_ntl_ulong)(posn+i))%NTL_BITS_PER_LONG] = -k+1;
         }
      }
   }
}

GF2XModulus& GF2XModulus::operator=(const GF2XModulus& F)
{
   if (this == &F) return *this;

   f=F.f; n=F.n; sn=F.sn; posn=F.posn; 
   k3=F.k3; k2=F.k2; k1=F.k1;
   size=F.size; 
   msk=F.msk; method=F.method; stab=F.stab; h0=F.h0; f0 = F.f0;
   tracevec=F.tracevec;

   if (method == GF2X_MOD_SPECIAL) {
      long i;
      if (!stab1) stab1 = NTL_NEW_OP _ntl_ulong[2*NTL_BITS_PER_LONG];
      if (!stab1) Error("GF2XModulus: out of memory");
      for (i = 0; i < 2*NTL_BITS_PER_LONG; i++)
         stab1[i] = F.stab1[i];
      if (!stab_cnt) stab_cnt = NTL_NEW_OP long[NTL_BITS_PER_LONG];
      if (!stab_cnt) Error("GF2XModulus: out of memory");
      for (i = 0; i < NTL_BITS_PER_LONG; i++)
         stab_cnt[i] = F.stab_cnt[i];
   }
   else if (method == GF2X_MOD_PLAIN) {
      long i;

      if (F.stab_cnt) {
         if (!stab_cnt) stab_cnt = NTL_NEW_OP long[NTL_BITS_PER_LONG];
         if (!stab_cnt) Error("GF2XModulus: out of memory");
         for (i = 0; i < NTL_BITS_PER_LONG; i++)
            stab_cnt[i] = F.stab_cnt[i];
      }

      if (F.stab_ptr) {
         if (!stab_ptr) stab_ptr = NTL_NEW_OP _ntl_ulong_ptr[NTL_BITS_PER_LONG];
         if (!stab_ptr) Error("GF2XModulus: out of memory");
      
         for (i = 0; i < NTL_BITS_PER_LONG; i++) {
            WordVector& st = stab[((_ntl_ulong)(posn+i))%NTL_BITS_PER_LONG].xrep;
            long k = st.length();
            stab_ptr[((_ntl_ulong)(posn+i))%NTL_BITS_PER_LONG] = &st[k-1];
            stab_cnt[((_ntl_ulong)(posn+i))%NTL_BITS_PER_LONG] = -k+1;
         }
      }
   }

   return *this;
}
   


GF2XModulus::~GF2XModulus() 
{ 
   delete [] stab_ptr; 
   delete [] stab_cnt; 
   delete [] stab1; 
}



GF2XModulus::GF2XModulus(const GF2X& ff)
{
   n = -1;
   method = GF2X_MOD_PLAIN;
   stab_ptr = 0;
   stab_cnt = 0;
   stab1 = 0;

   build(*this, ff);
}





void UseMulRem21(GF2X& r, const GF2X& a, const GF2XModulus& F)
{
   GF2XRegister(P1);
   GF2XRegister(P2);

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   add(P2, P2, P1);
   mul(P1, P2, F.f0);
   trunc(P1, P1, F.n);
   trunc(r, a, F.n);
   add(r, r, P1);
}

void UseMulDivRem21(GF2X& q, GF2X& r, const GF2X& a, const GF2XModulus& F)
{
   GF2XRegister(P1);
   GF2XRegister(P2);

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   add(P2, P2, P1);
   mul(P1, P2, F.f0);
   trunc(P1, P1, F.n);
   trunc(r, a, F.n);
   add(r, r, P1);
   q = P2;
}

void UseMulDiv21(GF2X& q, const GF2X& a, const GF2XModulus& F)
{
   GF2XRegister(P1);
   GF2XRegister(P2);

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   add(P2, P2, P1);
   q = P2;
}


void UseMulRemX1(GF2X& r, const GF2X& aa, const GF2XModulus& F)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);

   clear(buf);
   a = aa;

   long n = F.n;
   long a_len = deg(a) + 1;

   while (a_len > 0) {
      long old_buf_len = deg(buf) + 1;
      long amt = min(2*n-1-old_buf_len, a_len);

      LeftShift(buf, buf, amt);
      RightShift(tmp, a, a_len-amt);
      add(buf, buf, tmp);
      trunc(a, a, a_len-amt);

      UseMulRem21(buf, buf, F);
      a_len -= amt;
   }

   r = buf;
}
   

void UseMulDivRemX1(GF2X& q, GF2X& r, const GF2X& aa, const GF2XModulus& F)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);
   GF2XRegister(qq);
   GF2XRegister(qbuf);

   clear(buf);
   a = aa;
   clear(qq);

   long n = F.n;
   long a_len = deg(a) + 1;

   while (a_len > 0) {
      long old_buf_len = deg(buf) + 1;
      long amt = min(2*n-1-old_buf_len, a_len);

      LeftShift(buf, buf, amt);
      RightShift(tmp, a, a_len-amt);
      add(buf, buf, tmp);
      trunc(a, a, a_len-amt);

      UseMulDivRem21(qbuf, buf, buf, F);
      a_len -= amt;

      ShiftAdd(qq, qbuf, a_len);
   }

   r = buf;
   q = qq;
}


void UseMulDivX1(GF2X& q, const GF2X& aa, const GF2XModulus& F)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);
   GF2XRegister(qq);
   GF2XRegister(qbuf);
   
   clear(buf);
   a = aa;
   clear(qq);

   long n = F.n;
   long a_len = deg(a) + 1;

   while (a_len > 0) {
      long old_buf_len = deg(buf) + 1;
      long amt = min(2*n-1-old_buf_len, a_len);

      LeftShift(buf, buf, amt);
      RightShift(tmp, a, a_len-amt);
      add(buf, buf, tmp);
      trunc(a, a, a_len-amt);

      UseMulDivRem21(qbuf, buf, buf, F);
      a_len -= amt;