gf2x1.c

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

      ShiftAdd(qq, qbuf, a_len);
   }

   q = qq;
}

static
void TrinomReduce(GF2X& x, const GF2X& a, long n, long k)
{
   long wn = n / NTL_BITS_PER_LONG;
   long bn = n - wn*NTL_BITS_PER_LONG;

   long wdiff = (n-k)/NTL_BITS_PER_LONG;
   long bdiff = (n-k) - wdiff*NTL_BITS_PER_LONG;

   long m = a.xrep.length()-1;

   if (wn > m) {
      x = a;
      return;
   }

   GF2XRegister(r);

   r = a;

   _ntl_ulong *p = r.xrep.elts();

   _ntl_ulong *pp;


   _ntl_ulong w;

   if (bn == 0) {
      if (bdiff == 0) {
         // bn == 0 && bdiff == 0

         while (m >= wn) {
            w = p[m];
            p[m-wdiff] ^= w;
            p[m-wn] ^= w;
            m--;
         }
      }
      else {
         // bn == 0 && bdiff != 0

         while (m >= wn) {
            w = p[m];
            pp = &p[m-wdiff];
            *pp ^= (w >> bdiff);
            *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bdiff));
            p[m-wn] ^= w;
            m--;
         }
      }
   }
   else {
      if (bdiff == 0) {
         // bn != 0 && bdiff == 0

         while (m > wn) {
            w = p[m];
            p[m-wdiff] ^= w;
            pp = &p[m-wn];
            *pp ^= (w >> bn);
            *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bn));
            m--;
         }

         w = (p[m] >> bn) << bn;;

         p[m-wdiff] ^= w;
         p[0] ^= (w >> bn);

         p[m] &= ((1UL<<bn)-1UL); 
      }
      else {
         // bn != 0 && bdiff != 0

         while (m > wn) {
            w = p[m];
            pp = &p[m-wdiff];
            *pp ^= (w >> bdiff);;
            *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bdiff));
            pp = &p[m-wn];
            *pp ^= (w >> bn);
            *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bn));
            m--;
         }

         w = (p[m] >> bn) << bn;;

         p[m-wdiff] ^= (w >> bdiff);
         if (m-wdiff-1 >= 0) p[m-wdiff-1] ^= (w << (NTL_BITS_PER_LONG-bdiff));
         p[0] ^= (w >> bn);
         p[m] &= ((1UL<<bn)-1UL); 
      }
   }

   if (bn == 0)
      wn--;

   while (wn >= 0 && p[wn] == 0)
      wn--;

   r.xrep.QuickSetLength(wn+1);

   x = r;
}

static
void PentReduce(GF2X& x, const GF2X& a, long n, long k3, long k2, long k1)
{
   long wn = n / NTL_BITS_PER_LONG;
   long bn = n - wn*NTL_BITS_PER_LONG;

   long m = a.xrep.length()-1;

   if (wn > m) {
      x = a;
      return;
   }

   long wdiff1 = (n-k1)/NTL_BITS_PER_LONG;
   long bdiff1 = (n-k1) - wdiff1*NTL_BITS_PER_LONG;

   long wdiff2 = (n-k2)/NTL_BITS_PER_LONG;
   long bdiff2 = (n-k2) - wdiff2*NTL_BITS_PER_LONG;

   long wdiff3 = (n-k3)/NTL_BITS_PER_LONG;
   long bdiff3 = (n-k3) - wdiff3*NTL_BITS_PER_LONG;

   static GF2X r;
   r = a;

   _ntl_ulong *p = r.xrep.elts();

   _ntl_ulong *pp;

   _ntl_ulong w;

   while (m > wn) {
      w = p[m];

      if (bn == 0) 
         p[m-wn] ^= w;
      else {
         pp = &p[m-wn];
         *pp ^= (w >> bn);
         *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bn));
      }

      if (bdiff1 == 0) 
         p[m-wdiff1] ^= w;
      else {
         pp = &p[m-wdiff1];
         *pp ^= (w >> bdiff1);
         *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bdiff1));
      }

      if (bdiff2 == 0) 
         p[m-wdiff2] ^= w;
      else {
         pp = &p[m-wdiff2];
         *pp ^= (w >> bdiff2);
         *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bdiff2));
      }

      if (bdiff3 == 0) 
         p[m-wdiff3] ^= w;
      else {
         pp = &p[m-wdiff3];
         *pp ^= (w >> bdiff3);
         *(pp-1) ^= (w << (NTL_BITS_PER_LONG-bdiff3));
      }

      m--;
   }

   w = (p[m] >> bn) << bn;

   p[0] ^= (w >> bn); 

   if (bdiff1 == 0)
      p[m-wdiff1] ^= w;
   else {
      p[m-wdiff1] ^= (w >> bdiff1);
      if (m-wdiff1-1 >= 0) p[m-wdiff1-1] ^= (w << (NTL_BITS_PER_LONG-bdiff1));
   }

   if (bdiff2 == 0)
      p[m-wdiff2] ^= w;
   else {
      p[m-wdiff2] ^= (w >> bdiff2);
      if (m-wdiff2-1 >= 0) p[m-wdiff2-1] ^= (w << (NTL_BITS_PER_LONG-bdiff2));
   }

   if (bdiff3 == 0)
      p[m-wdiff3] ^= w;
   else {
      p[m-wdiff3] ^= (w >> bdiff3);
      if (m-wdiff3-1 >= 0) p[m-wdiff3-1] ^= (w << (NTL_BITS_PER_LONG-bdiff3));
   }

   if (bn != 0)
      p[m] &= ((1UL<<bn)-1UL);

   
   if (bn == 0)
      wn--;

   while (wn >= 0 && p[wn] == 0)
      wn--;

   r.xrep.QuickSetLength(wn+1);

   x = r;
}




static
void RightShiftAdd(GF2X& c, const GF2X& a, long n)
{
   if (n < 0) {
      Error("RightShiftAdd: negative shamt");
   }

   if (n == 0) {
      add(c, c, a);
      return;
   }

   long sa = a.xrep.length();
   long wn = n/NTL_BITS_PER_LONG;
   long bn = n - wn*NTL_BITS_PER_LONG;

   if (wn >= sa) {
      return;
   }

   long sc = c.xrep.length();
   long i;

   if (sa-wn > sc)
      c.xrep.SetLength(sa-wn);

   _ntl_ulong *cp = c.xrep.elts();
   const _ntl_ulong *ap = a.xrep.elts();

   for (i = sc; i < sa-wn; i++)
      cp[i] = 0;


   if (bn == 0) {
      for (i = 0; i < sa-wn; i++)
         cp[i] ^= ap[i+wn];
   }
   else {
      for (i = 0; i < sa-wn-1; i++)
         cp[i] ^= (ap[i+wn] >> bn) | (ap[i+wn+1] << (NTL_BITS_PER_LONG - bn));

      cp[sa-wn-1] ^= ap[sa-1] >> bn;
   }

   c.normalize();
}


static
void TriDiv21(GF2X& q, const GF2X& a, long n, long k)
{
   GF2XRegister(P1);

   RightShift(P1, a, n);
   if (k != 1) 
      RightShiftAdd(P1, P1, n-k);

   q = P1;
}

static 
void TriDivRem21(GF2X& q, GF2X& r, const GF2X& a, long n, long k)
{
   GF2XRegister(Q);
   TriDiv21(Q, a, n, k);
   TrinomReduce(r, a, n, k);
   q = Q;
}


static
void PentDiv21(GF2X& q, const GF2X& a, long n, long k3, long k2, long k1)
{
   if (deg(a) < n) {
      clear(q);
      return;
   }

   GF2XRegister(P1);
   GF2XRegister(P2);

   RightShift(P1, a, n);
   
   RightShift(P2, P1, n-k3);
   RightShiftAdd(P2, P1, n-k2);
   if (k1 != 1) {
      RightShiftAdd(P2, P1, n-k1);
   }

   add(P2, P2, P1);

   q = P2;
}

static 
void PentDivRem21(GF2X& q, GF2X& r, const GF2X& a, long n, 
                  long k3, long k2, long k1)
{
   GF2XRegister(Q);
   PentDiv21(Q, a, n, k3, k2, k1);
   PentReduce(r, a, n, k3, k2, k1);
   q = Q;
}

static
void TriDivRemX1(GF2X& q, GF2X& r, const GF2X& aa, long n, long k)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);
   GF2XRegister(qq);
   GF2XRegister(qbuf);

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

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

      TriDivRem21(qbuf, buf, buf, n, k);
      a_len -= amt;

      ShiftAdd(qq, qbuf, a_len);
   }

   r = buf;
   q = qq;
}


static
void TriDivX1(GF2X& q, const GF2X& aa, long n, long k)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);
   GF2XRegister(qq);
   GF2XRegister(qbuf);
   
   clear(buf);
   a = aa;
   clear(qq);

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

      TriDivRem21(qbuf, buf, buf, n, k);
      a_len -= amt;

      ShiftAdd(qq, qbuf, a_len);
   }

   q = qq;
}

static
void PentDivRemX1(GF2X& q, GF2X& r, const GF2X& aa, long n, 
                  long k3, long k2, long k1)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);
   GF2XRegister(qq);
   GF2XRegister(qbuf);

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

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

      PentDivRem21(qbuf, buf, buf, n, k3, k2, k1);
      a_len -= amt;

      ShiftAdd(qq, qbuf, a_len);
   }

   r = buf;
   q = qq;
}


static
void PentDivX1(GF2X& q, const GF2X& aa, long n, long k3, long k2, long k1)
{
   GF2XRegister(buf);
   GF2XRegister(tmp);
   GF2XRegister(a);
   GF2XRegister(qq);
   GF2XRegister(qbuf);
   
   clear(buf);
   a = aa;
   clear(qq);

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

      PentDivRem21(qbuf, buf, buf, n, k3, k2, k1);
      a_len -= amt;

      ShiftAdd(qq, qbuf, a_len);
   }

   q = qq;
}



void rem(GF2X& r, const GF2X& a, const GF2XModulus& F)
{
   long n = F.n;

   if (n < 0) Error("rem: uninitialized modulus");

   if (F.method == GF2X_MOD_TRI) {
      TrinomReduce(r, a, n, F.k3);
      return;
   }

   if (F.method == GF2X_MOD_PENT) {
      PentReduce(r, a, n, F.k3, F.k2, F.k1);
      return;
   }

   long da = deg(a);


   if (da < n) {
      r = a;
   }
   else if (F.method == GF2X_MOD_MUL) {
      if (da <= 2*(n-1)) 
         UseMulRem21(r, a, F);
      else
         UseMulRemX1(r, a, F);
   }
   else if (F.method == GF2X_MOD_SPECIAL) {
      long sa = a.xrep.length();
      long posa = da - NTL_BITS_PER_LONG*(sa-1);
   
      _ntl_ulong *ap;
      if (&r == &a)
         ap = r.xrep.elts();
      else {
         GF2X_rembuf = a.xrep;
         ap = GF2X_rembuf.elts();
      }
   
      _ntl_ulong *atop = &ap[sa-1];
      _ntl_ulong *stab_top;
   
      long i;
   
      while (da >= n) {
         if (atop[0] & (1UL << posa)) {
            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--;
         }
      }
   
      long sn = F.size;
      r.xrep.SetLength(sn);
      if (&r != &a) {
         _ntl_ulong *rp = r.xrep.elts();
         for (i = 0; i < sn; i++)
            rp[i] = ap[i];
      }
      r.xrep[sn-1] &= F.msk;
      r.normalize();
   }
   else {
      long sa = a.xrep.length();
      long posa = da - NTL_BITS_PER_LONG*(sa-1);
   
   
      _ntl_ulong *ap;
      if (&r == &a)
         ap = r.xrep.elts();
      else {
         GF2X_rembuf = a.xrep;
         ap = GF2X_rembuf.elts();
      }
   
      _ntl_ulong *atop = &ap[sa-1];
      _ntl_ulong *stab_top;
   
      long i;
   
      while (da >= n) {
         if (atop[0] & (1UL << posa)) {
            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--;
         }
      }
   
      long sn = F.size;
      r.xrep.SetLength(sn);
      if (&r != &a) {
         _ntl_ulong *rp = r.xrep.elts();
         for (i = 0; i < sn; i++)
            rp[i] = ap[i];
      }
      r.normalize();
   }
}

void DivRem(GF2X& q, GF2X& r, const GF2X& a, const GF2XModulus& F)
{
   long da = deg(a);
   long n = F.n;

   if (n < 0) Error("DivRem: uninitialized modulus");

   if (da < n) {
      r = a;
      clear(q);
   }
   else if (F.method == GF2X_MOD_TRI) {
      if (da <= 2*(n-1)) 
         TriDivRem21(q, r, a, F.n, F.k3);
      else
         TriDivRemX1(q, r, a, F.n, F.k3);
   }
   else if (F.method == GF2X_MOD_PENT) {
      if (da <= 2*(n-1)) 
         PentDivRem21(q, r, a, F.n, F.k3, F.k2, F.k1);
      else
         PentDivRemX1(q, r, a, F.n, F.k3, F.k2, F.k1);
   }
   else if (F.method == GF2X_MOD_MUL) {
      if (da <= 2*(n-1)) 
         UseMulDivRem21(q, r, a, F);
      else
         UseMulDivRemX1(q, r, a, F);
   }
   else if (F.method == GF2X_MOD_SPECIAL) {
      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;
      if (&r == &a)
         ap = r.xrep.elts();
      else {
         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.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--;
         }
      }
   
      long sn = F.size;
      r.xrep.SetLength(sn);
      if (&r != &a) {
         _ntl_ulong *rp = r.xrep.elts();
         for (i = 0; i < sn; i++)
            rp[i] = ap[i];
      }
      r.xrep[sn-1] &= F.msk;
      r.normalize();
   }
   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;
      if (&r == &a)
         ap = r.xrep.elts();
      else {
         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--;
         }
      }
   
      long sn = F.size;
      r.xrep.SetLength(sn);
      if (&r != &a) {
         _ntl_ulong *rp = r.xrep.elts();
         for (i = 0; i < sn; i++)
            rp[i] = ap[i];
      }
      r.normalize();
   }
}


void div(GF2X& q, const GF2X& a, const GF2XModulus& F)
{
   long da = deg(a);
   long n = F.n;

   if (n < 0) Error("div: uninitialized modulus");


   if (da < n) {
      clear(q);
   }
   else if (F.method == GF2X_MOD_TRI) {
      if (da <= 2*(n-1)) 
         TriDiv21(q, a, F.n, F.k3);
      else
         TriDivX1(q, a, F.n, F.k3);
   }
   else if (F.method == GF2X_MOD_PENT) {
      if (da <= 2*(n-1)) 
         PentDiv21(q, a, F.n, F.k3, F.k2, F.k1);
      else
         PentDivX1(q, a, F.n, F.k3, F.k2, F.k1);
   }
   else if (F.method == GF2X_MOD_MUL) {
      if (da <= 2*(n-1)) 
         UseMulDiv21(q, a, F);
      else
         UseMulDivX1(q, a, F);
   }
   else if (F.method == GF2X_MOD_SPECIAL) {
      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();