gf2x1.c

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

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

      long i;

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

      _ntl_ulong *qtop = &qp[sq-1];
   
      while (da >= n) {
         if (atop[0] & (1UL << posa)) {
            qtop[0] |= (1UL << posq);
            stab_top = &F.stab1[posa << 1];
            i = F.stab_cnt[posa];
            atop[i] ^= stab_top[0];
            atop[i+1] ^= stab_top[1];
         }
   
         da--;
         posa--;
         if (posa < 0) {
            posa = NTL_BITS_PER_LONG-1;
            atop--;
         }

         posq--;
         if (posq < 0) {
            posq = NTL_BITS_PER_LONG-1;
            qtop--;
         }
      }
   }
   else {
      long sa = a.xrep.length();
      long posa = da - NTL_BITS_PER_LONG*(sa-1);
   
      long dq = da - n;
      long sq = dq/NTL_BITS_PER_LONG + 1;
      long posq = dq - NTL_BITS_PER_LONG*(sq-1);
   
      _ntl_ulong *ap;
      GF2X_rembuf = a.xrep;
      ap = GF2X_rembuf.elts();
   
      _ntl_ulong *atop = &ap[sa-1];
      _ntl_ulong *stab_top;
   
      long i;

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

      _ntl_ulong *qtop = &qp[sq-1];
   
      while (da >= n) {
         if (atop[0] & (1UL << posa)) {
            qtop[0] |= (1UL << posq);
            stab_top = F.stab_ptr[posa];
            for (i = F.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--;
         }
      }
   }
}


void MulMod(GF2X& c, const GF2X& a, const GF2X& b, const GF2XModulus& F)
{
   if (F.n < 0) Error("MulMod: uninitialized modulus");

   GF2XRegister(t);
   mul(t, a, b);
   rem(c, t, F);
}


void SqrMod(GF2X& c, const GF2X& a, const GF2XModulus& F)
{
   if (F.n < 0) Error("SqrMod: uninitialized modulus");

   GF2XRegister(t);
   sqr(t, a);
   rem(c, t, F);
}


// we need these two versions to prevent a GF2XModulus
// from being constructed.


void MulMod(GF2X& c, const GF2X& a, const GF2X& b, const GF2X& f)
{
   GF2XRegister(t);
   mul(t, a, b);
   rem(c, t, f);
}

void SqrMod(GF2X& c, const GF2X& a, const GF2X& f)
{
   GF2XRegister(t);
   sqr(t, a);
   rem(c, t, f);
}


static
long OptWinSize(long n)
// finds k that minimizes n/(k+1) + 2^{k-1}

{
   long k;
   double v, v_new;


   v = n/2.0 + 1.0;
   k = 1;

   for (;;) {
      v_new = n/(double(k+2)) + double(1L << k);
      if (v_new >= v) break;
      v = v_new;
      k++;
   }

   return k;
}
      


void PowerMod(GF2X& h, const GF2X& g, const ZZ& e, const GF2XModulus& F)
// h = g^e mod f using "sliding window" algorithm
{
   if (deg(g) >= F.n) Error("PowerMod: bad args");

   if (e == 0) {
      set(h);
      return;
   }

   if (e == 1) {
      h = g;
      return;
   }

   if (e == -1) {
      InvMod(h, g, F);
      return;
   }

   if (e == 2) {
      SqrMod(h, g, F);
      return;
   }

   if (e == -2) {
      SqrMod(h, g, F);
      InvMod(h, h, F);
      return;
   }


   long n = NumBits(e);

   GF2X res;
   res.SetMaxLength(F.n);
   set(res);

   long i;

   if (n < 16) {
      // plain square-and-multiply algorithm

      for (i = n - 1; i >= 0; i--) {
         SqrMod(res, res, F);
         if (bit(e, i))
            MulMod(res, res, g, F);
      }

      if (e < 0) InvMod(res, res, F);

      h = res;
      return;
   }

   long k = OptWinSize(n);

   k = min(k, 9);

   vec_GF2X v;

   v.SetLength(1L << (k-1));

   v[0] = g;
 
   if (k > 1) {
      GF2X t;
      SqrMod(t, g, F);

      for (i = 1; i < (1L << (k-1)); i++)
         MulMod(v[i], v[i-1], t, F);
   }


   long val;
   long cnt;
   long m;

   val = 0;
   for (i = n-1; i >= 0; i--) {
      val = (val << 1) | bit(e, i); 
      if (val == 0)
         SqrMod(res, res, F);
      else if (val >= (1L << (k-1)) || i == 0) {
         cnt = 0;
         while ((val & 1) == 0) {
            val = val >> 1;
            cnt++;
         }

         m = val;
         while (m > 0) {
            SqrMod(res, res, F);
            m = m >> 1;
         }

         MulMod(res, res, v[val >> 1], F);

         while (cnt > 0) {
            SqrMod(res, res, F);
            cnt--;
         }

         val = 0;
      }
   }

   if (e < 0) InvMod(res, res, F);

   h = res;
}

   


void PowerXMod(GF2X& hh, const ZZ& e, const GF2XModulus& F)
{
   if (F.n < 0) Error("PowerXMod: uninitialized modulus");

   if (IsZero(e)) {
      set(hh);
      return;
   }

   long n = NumBits(e);
   long i;

   GF2X h;

   h.SetMaxLength(F.n+1);
   set(h);

   for (i = n - 1; i >= 0; i--) {
      SqrMod(h, h, F);
      if (bit(e, i)) {
         MulByX(h, h);
         if (coeff(h, F.n) != 0)
            add(h, h, F.f);
      }
   }

   if (e < 0) InvMod(h, h, F);

   hh = h;
}


      


void UseMulRem(GF2X& r, const GF2X& a, const GF2X& b)
{
   GF2XRegister(P1);
   GF2XRegister(P2);

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

   CopyReverse(P1, b, db);
   InvTrunc(P2, P1, da-db+1);
   CopyReverse(P1, P2, da-db);

   RightShift(P2, a, db);
   mul(P2, P1, P2);
   RightShift(P2, P2, da-db);
   mul(P1, P2, b);
   add(P1, P1, a);
   
   r = P1;
}

void UseMulDivRem(GF2X& q, GF2X& r, const GF2X& a, const GF2X& b)
{
   GF2XRegister(P1);
   GF2XRegister(P2);

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

   CopyReverse(P1, b, db);
   InvTrunc(P2, P1, da-db+1);
   CopyReverse(P1, P2, da-db);

   RightShift(P2, a, db);
   mul(P2, P1, P2);
   RightShift(P2, P2, da-db);
   mul(P1, P2, b);
   add(P1, P1, a);
   
   r = P1;
   q = P2;
}

void UseMulDiv(GF2X& q, const GF2X& a, const GF2X& b)
{
   GF2XRegister(P1);
   GF2XRegister(P2);

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

   CopyReverse(P1, b, db);
   InvTrunc(P2, P1, da-db+1);
   CopyReverse(P1, P2, da-db);

   RightShift(P2, a, db);
   mul(P2, P1, P2);
   RightShift(P2, P2, da-db);
   
   q = P2;
}


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


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



void GCD(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;
   }
}


#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];  \
         }  \
      }  \





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


void XGCD(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 InvMod(GF2X& c, const GF2X& a, const GF2X& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args");
   
   GF2XRegister(d);
   GF2XRegister(s);
   XXGCD(d, s, a, f);

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

   c = s;
}

long InvModStatus(GF2X& c, const GF2X& a, const GF2X& f)
{
   if (deg(a) >= deg(f) || deg(f) == 0) Error("InvMod: bad args");
   
   GF2XRegister(d);
   GF2XRegister(s);
   XXGCD(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 (n >= (1L << (NTL_BITS_PER_LONG-4)))
      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();