gf2ex.cpp

数值算法库for Windows

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

   GF2EX t1, t2;

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

   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);
   LeftShift(t2, t2, 1);

   add(x, t1, t2);
}


void UpdateMap(vec_GF2E& x, const vec_GF2E& a, 
         const GF2EXTransMultiplier& B, const GF2EXModulus& F)
{
   GF2EX xx;
   TransMulMod(xx, to_GF2EX(a), B, F);
   x = xx.rep;
}
   


static
void ProjectPowers(vec_GF2E& x, const GF2EX& a, long k, 
                   const GF2EXArgument& H, const GF2EXModulus& F)
{
   if (k < 0 || k >= (1L << (NTL_BITS_PER_LONG-4)) || deg(a) >= F.n) 
      Error("ProjectPowers: bad args");

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

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

   GF2EX s;
   s = a;

   x.SetLength(k);

   long i;

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

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

   if (k == 0) {
      x.SetLength(0);;
      return;
   }

   long m = SqrRoot(k);

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

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

void ProjectPowers(vec_GF2E& x, const vec_GF2E& a, long k,
                   const GF2EXArgument& H, const GF2EXModulus& F)
{
   ProjectPowers(x, to_GF2EX(a), k, H, F);
}

void ProjectPowers(vec_GF2E& x, const vec_GF2E& a, long k, 
                   const GF2EX& h, const GF2EXModulus& F)
{
   ProjectPowers(x, to_GF2EX(a), k, h, F);
}




void BerlekampMassey(GF2EX& h, const vec_GF2E& a, long m)
{
   GF2EX Lambda, Sigma, Temp;
   long L;
   GF2E Delta, Delta1, t1;
   long shamt;
   GF2X tt1, tt2;

   // cerr << "*** " << m << "\n";

   Lambda.SetMaxLength(m+1);
   Sigma.SetMaxLength(m+1);
   Temp.SetMaxLength(m+1);

   L = 0;
   set(Lambda);
   clear(Sigma);
   set(Delta);
   shamt = 0;

   long i, r, dl;

   for (r = 1; r <= 2*m; r++) {
      // cerr << r << "--";
      clear(tt1);
      dl = deg(Lambda);
      for (i = 0; i <= dl; i++) {
         mul(tt2, rep(Lambda.rep[i]), rep(a[r-i-1]));
         add(tt1, tt1, tt2);
      }

      conv(Delta1, tt1);

      if (IsZero(Delta1)) {
         shamt++;
         // cerr << "case 1: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else if (2*L < r) {
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         Sigma = Lambda;
         ShiftAdd(Lambda, Temp, shamt+1);
         shamt = 0;
         L = r-L;
         Delta = Delta1;
         // cerr << "case 2: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
      else {
         shamt++;
         div(t1, Delta1, Delta);
         mul(Temp, Sigma, t1);
         ShiftAdd(Lambda, Temp, shamt);
         // cerr << "case 3: " << deg(Lambda) << " " << deg(Sigma) << " " << shamt << "\n";
      }
   }

   // cerr << "finished: " << L << " " << deg(Lambda) << "\n"; 

   dl = deg(Lambda);
   h.rep.SetLength(L + 1);

   for (i = 0; i < L - dl; i++)
      clear(h.rep[i]);

   for (i = L - dl; i <= L; i++)
      h.rep[i] = Lambda.rep[L - i];
}


void MinPolySeq(GF2EX& h, const vec_GF2E& 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");

   BerlekampMassey(h, a, m);
}


void DoMinPolyMod(GF2EX& h, const GF2EX& g, const GF2EXModulus& F, long m, 
               const GF2EX& R)
{