gf2x1.cpp

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

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;

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