gf2ex.c

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

      PlainDiv(q, a, F.f);
      return;
   }

   long da = deg(a);
   long n = F.n;

   if (da <= 2*n-2) {
      UseMulDiv21(q, a, F);
      return;
   }

   GF2EX buf(INIT_SIZE, 2*n-1);
   GF2EX qbuf(INIT_SIZE, n-1);

   GF2EX qq;
   qq.rep.SetLength(da-n+1);

   long a_len = da+1;
   long q_hi = da-n+1;

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

      buf.rep.SetLength(old_buf_len+amt);

      long i;

      for (i = old_buf_len+amt-1; i >= amt; i--)
         buf.rep[i] = buf.rep[i-amt];

      for (i = amt-1; i >= 0; i--)
         buf.rep[i] = a.rep[a_len-amt+i];

      buf.normalize();

      a_len = a_len - amt;
      if (a_len > 0)
         UseMulDivRem21(qbuf, buf, buf, F);
      else
         UseMulDiv21(qbuf, buf, F);

      long dl = qbuf.rep.length();
      for(i = 0; i < dl; i++)
         qq.rep[a_len+i] = qbuf.rep[i];
      for(i = dl+a_len; i < q_hi; i++)
         clear(qq.rep[i]);
      q_hi = a_len;
   }

   qq.normalize();
   q = qq;
}




void MulMod(GF2EX& c, const GF2EX& a, const GF2EX& b, const GF2EXModulus& F)
{
   if (deg(a) >= F.n || deg(b) >= F.n) Error("MulMod: bad args");

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


void SqrMod(GF2EX& c, const GF2EX& a, const GF2EXModulus& F)
{
   if (deg(a) >= F.n) Error("MulMod: bad args");

   GF2EX 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(GF2EX& h, const GF2EX& g, const ZZ& e, const GF2EXModulus& 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);

   GF2EX 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, 5);

   vec_GF2EX v;

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

   v[0] = g;
 
   if (k > 1) {
      GF2EX 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(GF2EX& hh, const ZZ& e, const GF2EXModulus& F)
{
   if (F.n < 0) Error("PowerXMod: uninitialized modulus");

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

   long n = NumBits(e);
   long i;

   GF2EX h;

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

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

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

   hh = h;
}


      


void UseMulRem(GF2EX& r, const GF2EX& a, const GF2EX& b)
{
   GF2EX P1;
   GF2EX 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(GF2EX& q, GF2EX& r, const GF2EX& a, const GF2EX& b)
{
   GF2EX P1;
   GF2EX 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(GF2EX& q, const GF2EX& a, const GF2EX& b)
{
   GF2EX P1;
   GF2EX 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;
}



void DivRem(GF2EX& q, GF2EX& r, const GF2EX& a, const GF2EX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sb < GF2E::DivCross() || sa-sb < GF2E::DivCross())
      PlainDivRem(q, r, a, b);
   else if (sa < 4*sb)
      UseMulDivRem(q, r, a, b);
   else {
      GF2EXModulus B;
      build(B, b);
      DivRem(q, r, a, B);
   }
}

void div(GF2EX& q, const GF2EX& a, const GF2EX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sb < GF2E::DivCross() || sa-sb < GF2E::DivCross())
      PlainDiv(q, a, b);
   else if (sa < 4*sb)
      UseMulDiv(q, a, b);
   else {
      GF2EXModulus B;
      build(B, b);
      div(q, a, B);
   }
}

void div(GF2EX& q, const GF2EX& a, const GF2E& b)
{
   GF2E t;
   inv(t, b);
   mul(q, a, t);
}

void div(GF2EX& q, const GF2EX& a, GF2 b)
{
   if (b == 0)
      Error("div: division by zero");

   q = a;
}

void div(GF2EX& q, const GF2EX& a, long b)
{
   if ((b & 1) == 0)
      Error("div: division by zero");

   q = a;
}
   


void rem(GF2EX& r, const GF2EX& a, const GF2EX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

   if (sb < GF2E::DivCross() || sa-sb < GF2E::DivCross())
      PlainRem(r, a, b);
   else if (sa < 4*sb)
      UseMulRem(r, a, b);
   else {
      GF2EXModulus B;
      build(B, b);
      rem(r, a, B);
   }
}


void diff(GF2EX& x, const GF2EX& a)
{
   long n = deg(a);
   long i;

   if (n <= 0) {
      clear(x);
      return;
   }

   if (&x != &a)
      x.rep.SetLength(n);

   for (i = 0; i <= n-1; i++) {
      if ((i+1)&1)
         x.rep[i] = a.rep[i+1];
      else
         clear(x.rep[i]);
   }

   if (&x == &a)
      x.rep.SetLength(n);

   x.normalize();
}


void RightShift(GF2EX& x, const GF2EX& a, long n)
{
   if (n < 0) {
      if (n < -NTL_MAX_LONG) Error("overflow in RightShift");
      LeftShift(x, a, -n);
      return;
   }

   long da = deg(a);
   long i;
 
   if (da < n) {
      clear(x);
      return;
   }

   if (&x != &a)
      x.rep.SetLength(da-n+1);

   for (i = 0; i <= da-n; i++)
      x.rep[i] = a.rep[i+n];

   if (&x == &a)
      x.rep.SetLength(da-n+1);

   x.normalize();
}

void LeftShift(GF2EX& x, const GF2EX& a, long n)
{
   if (n < 0) {
      if (n < -NTL_MAX_LONG) Error("overflow in LeftShift");
      RightShift(x, a, -n);
      return;
   }

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

   if (IsZero(a)) {
      clear(x);
      return;
   }

   long m = a.rep.length();

   x.rep.SetLength(m+n);

   long i;
   for (i = m-1; i >= 0; i--)
      x.rep[i+n] = a.rep[i];

   for (i = 0; i < n; i++)
      clear(x.rep[i]);
}


void ShiftAdd(GF2EX& U, const GF2EX& V, long n)
// assumes input does not alias output
{
   if (IsZero(V))
      return;

   long du = deg(U);
   long dv = deg(V);

   long d = max(du, n+dv);

   U.rep.SetLength(d+1);
   long i;

   for (i = du+1; i <= d; i++)
      clear(U.rep[i]);

   for (i = 0; i <= dv; i++)
      add(U.rep[i+n], U.rep[i+n], V.rep[i]);

   U.normalize();
}

NTL_vector_impl(GF2EX,vec_GF2EX)

NTL_eq_vector_impl(GF2EX,vec_GF2EX)

NTL_io_vector_impl(GF2EX,vec_GF2EX)




void IterBuild(GF2E* a, long n)
{
   long i, k;
   GF2E b, t;

   if (n <= 0) return;

   for (k = 1; k <= n-1; k++) {
      b = a[k];
      add(a[k], b, a[k-1]);
      for (i = k-1; i >= 1; i--) {
         mul(t, a[i], b);
         add(a[i], t, a[i-1]);
      }
      mul(a[0], a[0], b);
   }
} 



void BuildFromRoots(GF2EX& x, const vec_GF2E& a)
{
   long n = a.length();

   if (n == 0) {
      set(x);
      return;
   }

   x.rep.SetMaxLength(n+1);
   x.rep = a;
   IterBuild(&x.rep[0], n);
   x.rep.SetLength(n+1);
   SetCoeff(x, n);
}



void eval(GF2E& b, const GF2EX& f, const GF2E& a)
// does a Horner evaluation
{
   GF2E acc;
   long i;

   clear(acc);
   for (i = deg(f); i >= 0; i--) {
      mul(acc, acc, a);
      add(acc, acc, f.rep[i]);
   }

   b = acc;
}



void eval(vec_GF2E& b, const GF2EX& f, const vec_GF2E& a)
// naive algorithm:  repeats Horner
{
   if (&b == &f.rep) {
      vec_GF2E bb;
      eval(bb, f, a);
      b = bb;
      return;
   }

   long m = a.length();
   b.SetLength(m);
   long i;
   for (i = 0; i < m; i++) 
      eval(b[i], f, a[i]);
}




void interpolate(GF2EX& f, const vec_GF2E& a, const vec_GF2E& b)
{
   long m = a.length();
   if (b.length() != m) Error("interpolate: vector length mismatch");

   if (m == 0) {
      clear(f);
      return;
   }

   vec_GF2E prod;
   prod = a;

   GF2E t1, t2;

   long k, i;

   vec_GF2E res;
   res.SetLength(m);

   for (k = 0; k < m; k++) {

      const GF2E& aa = a[k];

      set(t1);
      for (i = k-1; i >= 0; i--) {
         mul(t1, t1, aa);
         add(t1, t1, prod[i]);
      }

      clear(t2);
      for (i = k-1; i >= 0; i--) {
         mul(t2, t2, aa);
         add(t2, t2, res[i]);
      }


      inv(t1, t1);
      sub(t2, b[k], t2);
      mul(t1, t1, t2);

      for (i = 0; i < k; i++) {
         mul(t2, prod[i], t1);
         add(res[i], res[i], t2);
      }

      res[k] = t1;

      if (k < m-1) {
         if (k == 0)
            negate(prod[0], prod[0]);
         else {
            negate(t1, a[k]);
            add(prod[k], t1, prod[k-1]);
            for (i = k-1; i >= 1; i--) {
               mul(t2, prod[i], t1);
               add(prod[i], t2, prod[i-1]);
            }
            mul(prod[0], prod[0], t1);
         }
      }
   }

   while (m > 0 && IsZero(res[m-1])) m--;
   res.SetLength(m);
   f.rep = res;
}

   
void InnerProduct(GF2EX& x, const vec_GF2E& v, long low, long high, 
                   const vec_GF2EX& H, long n, GF2XVec& t)
{
   GF2X s;
   long i, j;

   for (j = 0; j < n; j++)
      clear(t[j]);

   high = min(high, v.length()-1);
   for (i = low; i <= high; i++) {
      const vec_GF2E& h = H[i-low].rep;
      long m = h.length();
      const GF2X& w = rep(v[i]);

      for (j = 0; j < m; j++) {
         mul(s, w, rep(h[j]));
         add(t[j], t[j], s);
      }
   }

   x.rep.SetLength(n);
   for (j = 0; j < n; j++)
      conv(x.rep[j], t[j]);
   x.normalize();
}


void CompMod(GF2EX& x, const GF2EX& g, const GF2EXArgument& A, 
             const GF2EXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   GF2EX s, t;
   GF2XVec scratch(F.n, 2*GF2E::WordLength());

   long m = A.H.length() - 1;
   long l = ((g.rep.length()+m-1)/m) - 1;

   const GF2EX& M = A.H[m];

   InnerProduct(t, g.rep, l*m, l*m + m - 1, A.H, F.n, scratch);
   for (long i = l-1; i >= 0; i--) {
      InnerProduct(s, g.rep, i*m, i*m + m - 1, A.H, F.n, scratch);
      MulMod(t, t, M, F);
      add(t, t, s);
   }

   x = t;
}


void build(GF2EXArgument& A, const GF2EX& h, const GF2EXModulus& F, long m)
{
   long i;

   if (m <= 0 || deg(h) >= F.n)
      Error("build GF2EXArgument: bad args");

   if (m > F.n) m = F.n;

   if (GF2EXArgBound > 0) {
      double sz = GF2E::WordLength()+4;
      sz = sz*(sizeof (_ntl_ulong));
      sz = sz*F.n;
      sz = sz + NTL_VECTOR_HEADER_SIZE + sizeof(vec_GF2E);
      sz = sz/1024;
      m = min(m, long(GF2EXArgBound/sz));
      m = max(m, 1);
   }

   A.H.SetLength(m+1);

   set(A.H[0]);
   A.H[1] = h;
   for (i = 2; i <= m; i++) 
      MulMod(A.H[i], A.H[i-1], h, F);
}




long GF2EXArgBound = 0;


void CompMod(GF2EX& x, const GF2EX& g, const GF2EX& h, const GF2EXModulus& F)
   // x = g(h) mod f
{
   long m = SqrRoot(g.rep.length());

   if (m == 0) {
      clear(x);
      return;
   }

   GF2EXArgument A;

   build(A, h, F, m);

   CompMod(x, g, A, F);
}




void Comp2Mod(GF2EX& x1, GF2EX& x2, const GF2EX& g1, const GF2EX& g2,
              const GF2EX& h, const GF2EXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length());

   if (m == 0) {
      clear(x1);
      clear(x2);
      return;
   }

   GF2EXArgument A;

   build(A, h, F, m);

   GF2EX xx1, xx2;

   CompMod(xx1, g1, A, F);
   CompMod(xx2, g2, A, F);

   x1 = xx1;
   x2 = xx2;
}

void Comp3Mod(GF2EX& x1, GF2EX& x2, GF2EX& x3, 
              const GF2EX& g1, const GF2EX& g2, const GF2EX& g3,
              const GF2EX& h, const GF2EXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());

   if (m == 0) {
      clear(x1);
      clear(x2);
      clear(x3);
      return;
   }

   GF2EXArgument A;

   build(A, h, F, m);

   GF2EX xx1, xx2, xx3;

   CompMod(xx1, g1, A, F);
   CompMod(xx2, g2, A, F);
   CompMod(xx3, g3, A, F);

   x1 = xx1;
   x2 = xx2;
   x3 = xx3;
}





void build(GF2EXTransMultiplier& B, const GF2EX& b, const GF2EXModulus& F)
{
   long db = deg(b);

   if (db >= F.n) Error("build TransMultiplier: bad args");

   GF2EX t;

   LeftShift(t, b, F.n-1);
   div(t, t, F);

   // we optimize for low degree b

   long d;

   d = deg(t);
   if (d < 0)
      B.shamt_fbi = 0;
   else
      B.shamt_fbi = F.n-2 - d; 

   CopyReverse(B.fbi, t, d);

   // The following code optimizes the case when 
   // f = X^n + low degree poly

   trunc(t, F.f, F.n);
   d = deg(t);
   if (d < 0)
      B.shamt = 0;
   else
      B.shamt = d;

   CopyReverse(B.f0, t, d);

   if (db < 0)
      B.shamt_b = 0;
   else
      B.shamt_b = db;

   CopyReverse(B.b, b, db);
}

void TransMulMod(GF2EX& x, const GF2EX& a, const GF2EXTransMultiplier& B,
               const GF2EXModulus& F)
{
   if (deg(a) >= F.n) Error("TransMulMod: bad args");