gf2x1.c

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

   long i;
   for (i = 0; i < wmin; i++)
      xp[i] = ap[i];

   if (wa < wx) {
      for (i = wa; i < wx; i++)
         xp[i] = 0;
   }
   else {
      long p = n % NTL_BITS_PER_LONG;
      if (p != 0)
         xp[wx-1] &= (1UL << p) - 1UL;
   }
}


void add(GF2X& c, const GF2X& a, long b)
{
   c = a;
   if (b & 1) {
      long n = c.xrep.length();
      if (n == 0) 
         set(c);
      else {
         c.xrep[0] ^= 1;
         if (n == 1 && !c.xrep[0]) c.xrep.SetLength(0);
      }
   }
}

void add(GF2X& c, const GF2X& a, GF2 b)
{
   add(c, a, rep(b));
}


void MulTrunc(GF2X& c, const GF2X& a, const GF2X& b, long n)
{
   GF2XRegister(t);

   mul(t, a, b);
   trunc(c, t, n);
}

void SqrTrunc(GF2X& c, const GF2X& a, long n)
{
   GF2XRegister(t);

   sqr(t, a);
   trunc(c, t, n);
}


long divide(GF2X& q, const GF2X& a, const GF2X& b)
{
   if (IsZero(b)) {
      if (IsZero(a)) {
         clear(q);
         return 1;
      }
      else
         return 0;
   }

   GF2XRegister(lq);
   GF2XRegister(r);

   DivRem(lq, r, a, b);
   if (!IsZero(r)) return 0;
   q = lq;
   return 1;
}

long divide(const GF2X& a, const GF2X& b)
{
   if (IsZero(b)) return IsZero(a);
   GF2XRegister(r);
   rem(r, a, b);
   if (!IsZero(r)) return 0;
   return 1;
}



/*** modular composition routines and data structures ***/


void InnerProduct(GF2X& x, const GF2X& v, long dv, long low, long high, 
                   const vec_GF2X& H, long n, WordVector& t)
{
   long i, j;

   _ntl_ulong *tp = t.elts();

   for (i = 0; i < n; i++)
      tp[i] = 0;


   long w_low = low/NTL_BITS_PER_LONG;
   long b_low = low - w_low*NTL_BITS_PER_LONG;

   
   const _ntl_ulong *vp = &v.xrep[w_low];
   _ntl_ulong msk = 1UL << b_low;
   _ntl_ulong vv = *vp;

   high = min(high, dv);

   i = low;
   for (;;) {
      if (vv & msk) {
         const WordVector& h = H[i-low].xrep;
         long m = h.length();
         const _ntl_ulong *hp = h.elts();
         for (j = 0; j < m; j++)
            tp[j] ^= hp[j];
      }

      i++;
      if (i > high) break;

      msk = msk << 1;
      if (!msk) {
         msk = 1UL;
         vp++;
         vv = *vp;
      }
   }

   x.xrep = t;
   x.normalize();
}


void CompMod(GF2X& x, const GF2X& g, const GF2XArgument& A, const GF2XModulus& F)
{
   long dg = deg(g);
   if (dg <= 0) {
      x = g;
      return;
   }

   GF2X s, t;
   WordVector scratch(INIT_SIZE, F.size);

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

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

   x = t;
}

void build(GF2XArgument& A, const GF2X& h, const GF2XModulus& F, long m)
{
   if (m <= 0 || deg(h) >= F.n) Error("build GF2XArgument: bad args");

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

   long i;

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


void CompMod(GF2X& x, const GF2X& g, const GF2X& h, const GF2XModulus& F)
   // x = g(h) mod f
{
   long m = SqrRoot(deg(g)+1);

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

   GF2XArgument A;

   build(A, h, F, m);

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




void Comp2Mod(GF2X& x1, GF2X& x2, const GF2X& g1, const GF2X& g2,
              const GF2X& h, const GF2XModulus& F)

{
   long m = SqrRoot(deg(g1) + deg(g2) + 2);

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

   GF2XArgument A;

   build(A, h, F, m);

   GF2X xx1, xx2;

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

   x1 = xx1;
   x2 = xx2;
}

void Comp3Mod(GF2X& x1, GF2X& x2, GF2X& x3, 
              const GF2X& g1, const GF2X& g2, const GF2X& g3,
              const GF2X& h, const GF2XModulus& F)

{
   long m = SqrRoot(deg(g1) + deg(g2) + deg(g3) + 3);

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

   GF2XArgument A;

   build(A, h, F, m);

   GF2X 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(GF2XTransMultiplier& B, const GF2X& b, const GF2XModulus& F)
{
   long db = deg(b);

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

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

   if (F.method != GF2X_MOD_TRI && F.method != GF2X_MOD_PENT) {
   
      // 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(GF2X& x, const GF2X& a, const GF2XTransMultiplier& B,
               const GF2XModulus& F)
{
   if (deg(a) >= F.n) Error("TransMulMod: bad args");

   GF2XRegister(t1);
   GF2XRegister(t2);
   GF2XRegister(t3);

   mul(t1, a, B.b);
   RightShift(t1, t1, B.shamt_b);

   if (F.method == GF2X_MOD_TRI) {
      RightShift(t2, a, F.k3);
      add(t2, t2, a);
   }
   else if (F.method == GF2X_MOD_PENT) {
      RightShift(t2, a, F.k3);
      RightShift(t3, a, F.k2);
      add(t2, t2, t3);
      RightShift(t3, a, F.k1);
      add(t2, t2, t3);
      add(t2, t2, a);
   }
   else {
      mul(t2, a, B.f0);
      RightShift(t2, t2, B.shamt);
   }

   trunc(t2, t2, F.n-1);

   mul(t2, t2, B.fbi);
   if (B.shamt_fbi > 0) LeftShift(t2, t2, B.shamt_fbi);
   trunc(t2, t2, F.n-1);
   MulByX(t2, t2);

   add(x, t1, t2);
}

void UpdateMap(vec_GF2& x, const vec_GF2& a, const GF2XTransMultiplier& B,
       const GF2XModulus& F)
{
   GF2XRegister(xx);
   GF2XRegister(aa);
   conv(aa, a);
   TransMulMod(xx, aa, B, F);
   conv(x, xx);
}
   

void ProjectPowers(GF2X& x, const GF2X& a, long k, const GF2XArgument& H,
                   const GF2XModulus& F)
{
   long n = F.n;

   if (deg(a) >= n || k < 0 || k >= (1L << (NTL_BITS_PER_LONG-4))) 
      Error("ProjectPowers: bad args");

   long m = H.H.length()-1;
   long l = (k+m-1)/m - 1;

   GF2XTransMultiplier M;
   build(M, H.H[m], F);

   GF2X s;
   s = a;

   x.SetMaxLength(k);
   clear(x);

   long i;

   for (i = 0; i <= l; i++) {
      long m1 = min(m, k-i*m);
      for (long j = 0; j < m1; j++)
         SetCoeff(x, i*m+j, InnerProduct(H.H[j].xrep, s.xrep));
      if (i < l)
         TransMulMod(s, s, M, F);
   }
}


void ProjectPowers(vec_GF2& x, const vec_GF2& a, long k, 
                   const GF2XArgument& H, const GF2XModulus& F)
{
   GF2X xx;
   ProjectPowers(xx, to_GF2X(a), k, H, F);
   VectorCopy(x, xx, k);
}


void ProjectPowers(GF2X& x, const GF2X& a, long k, const GF2X& h, 
                   const GF2XModulus& F)
{
   if (deg(a) >= F.n || k < 0) Error("ProjectPowers: bad args");

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

   long m = SqrRoot(k);

   GF2XArgument H;
   build(H, h, F, m);

   ProjectPowers(x, a, k, H, F);
}

void ProjectPowers(vec_GF2& x, const vec_GF2& a, long k, const GF2X& H,
                   const GF2XModulus& F)
{
   GF2X xx;
   ProjectPowers(xx, to_GF2X(a), k, H, F);
   VectorCopy(x, xx, k);
}


void MinPolyInternal(GF2X& h, const GF2X& x, long m)
{
   GF2X a, b, r, s;
   GF2X a_in, b_in;

   if (IsZero(x)) {
      set(h);
      return;
   }

   clear(a_in);
   SetCoeff(a_in, 2*m);

   CopyReverse(b_in, x, 2*m-1);
      
   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;
   }

   sp[0] = 1;

   long sr = 0;
   long ss = 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;

         swap(rp, sp);
         swap(sr, ss);
      }

      // da >= db

      if (db < m) break;

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

   a.normalize();
   b.normalize();
   r.normalize();
   s.normalize();

   if (!parity) {
      h = s;
   }
   else {
      h = r;
   }
}


void DoMinPolyMod(GF2X& h, const GF2X& g, const GF2XModulus& F, long m, 
               const GF2X& R)
{
   GF2X x;

   ProjectPowers(x, R, 2*m, g, F);
   MinPolyInternal(h, x, m);
}

void MinPolySeq(GF2X& h, const vec_GF2& a, long m)
{
   if (m < 0 || m >= (1L << (NTL_BITS_PER_LONG-4))) Error("MinPoly: bad args");
   if (a.length() < 2*m) Error("MinPoly: sequence too short");
   GF2X x;
   x.xrep = a.rep;
   x.normalize();
   MinPolyInternal(h, x, m);
}

void ProbMinPolyMod(GF2X& h, const GF2X& g, const GF2XModulus& F, long m)
{
   long n = F.n;
   if (m < 1 || m > n) Error("ProbMinPoly: bad args");

   GF2X R;
   random(R, n);

   DoMinPolyMod(h, g, F, m, R);
}

void ProbMinPolyMod(GF2X& h, const GF2X& g, const GF2XModulus& F)
{
   ProbMinPolyMod(h, g, F, F.n);
}

void MinPolyMod(GF2X& hh, const GF2X& g, const GF2XModulus& F, long m)
{
   GF2X h, h1;
   long n = F.n;
   if (m < 1 || m > n) Error("MinPoly: bad args");

   /* probabilistically compute min-poly */

   ProbMinPolyMod(h, g, F, m);
   if (deg(h) == m) { hh = h; return; }
   CompMod(h1, h, g, F);
   if (IsZero(h1)) { hh = h; return; }

   /* not completely successful...must iterate */


   GF2X h2, h3;
   GF2X R;
   GF2XTransMultiplier H1;
   

   for (;;) {
      random(R, n);
      build(H1, h1, F);
      TransMulMod(R, R, H1, F);
      DoMinPolyMod(h2, g, F, m-deg(h), R);

      mul(h, h, h2);
      if (deg(h) == m) { hh = h; return; }
      CompMod(h3, h2, g, F);
      MulMod(h1, h3, h1, F);
      if (IsZero(h1)) { hh = h; return; }
   }
}

void IrredPolyMod(GF2X& h, const GF2X& g, const GF2XModulus& F, long m)
{
   if (m < 1 || m > F.n) Error("IrredPoly: bad args");

   GF2X R;
   set(R);

   DoMinPolyMod(h, g, F, m, R);
}



void IrredPolyMod(GF2X& h, const GF2X& g, const GF2XModulus& F)
{
   IrredPolyMod(h, g, F, F.n);
}



void MinPolyMod(GF2X& hh, const GF2X& g, const GF2XModulus& F)
{
   MinPolyMod(hh, g, F, F.n);
}



void MulByXMod(GF2X& c, const GF2X& a, const GF2XModulus& F)
{
   long da = deg(a);
   long df = deg(F);
   if (da >= df) Error("MulByXMod: bad args"); 

   MulByX(c, a);

   if (da >= 0 && da == df-1)
      add(c, c, F);
}

static
void MulByXModAux(GF2X& c, const GF2X& a, const GF2X& f)
{
   long da = deg(a);
   long df = deg(f);
   if (da >= df) Error("MulByXMod: bad args"); 

   MulByX(c, a);

   if (da >= 0 && da == df-1)
      add(c, c, f);
}

void MulByXMod(GF2X& h, const GF2X& a, const GF2X& f)
{
   if (&h == &f) {
      GF2X hh;
      MulByXModAux(hh, a, f);
      h = hh;
   }
   else
      MulByXModAux(h, a, f);
}




void power(GF2X& x, const GF2X& a, long e)
{
   if (e < 0) {
      Error("power: negative exponent");
   }

   if (e == 0) {
      x = 1;
      return;
   }

   if (a == 0 || a == 1) {
      x = a;
      return;
   }

   long da = deg(a);

   if (da > (NTL_MAX_LONG-1)/e)
      Error("overflow in power");

   GF2X res;
   res.SetMaxLength(da*e + 1);
   res = 1;
   
   long k = NumBits(e);
   long i;

   for (i = k - 1; i >= 0; i--) {
      sqr(res, res);
      if (bit(e, i))
         mul(res, res, a);
   }

   x = res;
}


static
void FastTraceVec(vec_GF2& S, const GF2XModulus& f)
{
   long n = deg(f);

   if (n <= 0) Error("TraceVec: bad args");

   GF2X x = reverse(-LeftShift(reverse(diff(reverse(f)), n-1), n-1)/f, n-1);

   VectorCopy(S, x, n);
   S.put(0, to_GF2(n));
}

static
void PlainTraceVec(vec_GF2& S, const GF2X& f)
{
   long n = deg(f);

   if (n <= 0) 
      Error("TraceVec: bad args");

   if (n == 0) {
      S.SetLength(0);
      return;
   }

   GF2X x = reverse(-LeftShift(reverse(diff(reverse(f)), n-1), n-1)/f, n-1);

   VectorCopy(S, x, n); 
   S.put(0, to_GF2(n));
}


void TraceVec(vec_GF2& S, const GF2X& f)
{
   PlainTraceVec(S, f);
}

static
void ComputeTraceVec(const GF2XModulus& F)
{
   vec_GF2& S = *((vec_GF2 *) &F.tracevec);

   if (S.length() > 0)
      return;

   if (F.method == GF2X_MOD_PLAIN) {
      PlainTraceVec(S, F.f);
   }
   else {
      FastTraceVec(S, F);
   }
}

void TraceMod(GF2& x, const GF2X& a, const GF2XModulus& F)
{
   long n = F.n;

   if (deg(a) >= n)
      Error("trace: bad args");

   if (F.tracevec.length() == 0) 
      ComputeTraceVec(F);

   project(x, F.tracevec, a);
}

void TraceMod(GF2& x, const GF2X& a, const GF2X& f)
{
   if (deg(a) >= deg(f) || deg(f) <= 0)
      Error("trace: bad args");

   project(x, TraceVec(f), a);
}


NTL_END_IMPL