gf2x1.cpp

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

      _ntl_ulong *atop = &ap[sa-1];
      _ntl_ulong *stab_top;
   
      long i;
   
      while (1) {
         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--;
         if (da < n) break;

         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 (1) {
         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--;
         if (da < n) break;

         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 (1) {
         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--;
         if (da < n) break;

         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 (1) {
         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--;
         if (da < n) break;

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

   
      _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 (1) {
         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--;
         if (da < n) break;

         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 (1) {
         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--;
         if (da < n) break;

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