gf2ex.c

数值算法库for Unix

   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)
{
   vec_GF2E x;

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

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

   GF2EX R;
   random(R, n);

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

void ProbMinPolyMod(GF2EX& h, const GF2EX& g, const GF2EXModulus& F)
{
   ProbMinPolyMod(h, g, F, F.n);
}

void MinPolyMod(GF2EX& hh, const GF2EX& g, const GF2EXModulus& F, long m)
{
   GF2EX 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 */


   GF2EX h2, h3;
   GF2EX R;
   GF2EXTransMultiplier 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(GF2EX& h, const GF2EX& g, const GF2EXModulus& F, long m)
{
   if (m < 1 || m > F.n) Error("IrredPoly: bad args");

   GF2EX R;
   set(R);

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



void IrredPolyMod(GF2EX& h, const GF2EX& g, const GF2EXModulus& F)
{
   IrredPolyMod(h, g, F, F.n);
}



void MinPolyMod(GF2EX& hh, const GF2EX& g, const GF2EXModulus& F)
{
   MinPolyMod(hh, g, F, F.n);
}


void MakeMonic(GF2EX& x)
{
   if (IsZero(x))
      return;

   if (IsOne(LeadCoeff(x)))
      return;

   GF2E t;

   inv(t, LeadCoeff(x));
   mul(x, x, t);
}


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

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

long divide(const GF2EX& a, const GF2EX& b)
{
   if (IsZero(b)) return IsZero(a);
   GF2EX lq, r;
   DivRem(lq, r, a, b);
   if (!IsZero(r)) return 0; 
   return 1;
}


long operator==(const GF2EX& a, long b)
{
   if (b & 1)
      return IsOne(a);
   else
      return IsZero(a);
}


long operator==(const GF2EX& a, GF2 b)
{
   if (b == 1)
      return IsOne(a);
   else
      return IsZero(a);
}

long operator==(const GF2EX& a, const GF2E& b)
{
   if (IsZero(b))
      return IsZero(a);

   if (deg(a) != 0)
      return 0;

   return a.rep[0] == b;
}



void power(GF2EX& x, const GF2EX& 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 == 0) {
      x = power(ConstTerm(a), e);
      return;
   }


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

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

void reverse(GF2EX& x, const GF2EX& a, long hi)
{
   if (hi < -1) Error("reverse: bad args");
   if (hi >= (1L << (NTL_BITS_PER_LONG-4))) Error("overflow in reverse");

   if (&x == &a) {
      GF2EX tmp;
      CopyReverse(tmp, a, hi);
      x = tmp;
   }
   else
      CopyReverse(x, a, hi);
}


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

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

   S.SetLength(n);
   S[0] = n;

   long i;
   for (i = 1; i < n; i++)
      S[i] = coeff(x, i);
}


void PlainTraceVec(vec_GF2E& S, const GF2EX& ff)
{
   if (deg(ff) <= 0)
      Error("TraceVec: bad args");

   GF2EX f;
   f = ff;

   MakeMonic(f);

   long n = deg(f);

   S.SetLength(n);

   if (n == 0)
      return;

   long k, i;
   GF2X acc, t;
   GF2E t1;

   S[0] = n;

   for (k = 1; k < n; k++) {
      mul(acc, rep(f.rep[n-k]), k);

      for (i = 1; i < k; i++) {
         mul(t, rep(f.rep[n-i]), rep(S[k-i]));
         add(acc, acc, t);
      }

      conv(t1, acc);
      negate(S[k], t1);
   }
}

void TraceVec(vec_GF2E& S, const GF2EX& f)
{
   if (deg(f) < GF2E::DivCross())
      PlainTraceVec(S, f);
   else
      FastTraceVec(S, f);
}

static
void ComputeTraceVec(const GF2EXModulus& F)
{
   vec_GF2E& S = *((vec_GF2E *) &F.tracevec);

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

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

void TraceMod(GF2E& x, const GF2EX& a, const GF2EXModulus& F)
{
   long n = F.n;

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

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

   InnerProduct(x, a.rep, F.tracevec);
}

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

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


void PlainResultant(GF2E& rres, const GF2EX& a, const GF2EX& b)
{
   GF2E res;
 
   if (IsZero(a) || IsZero(b))
      clear(res);
   else if (deg(a) == 0 && deg(b) == 0) 
      set(res);
   else {
      long d0, d1, d2;
      GF2E lc;
      set(res);

      long n = max(deg(a),deg(b)) + 1;
      GF2EX u(INIT_SIZE, n), v(INIT_SIZE, n);
      GF2XVec tmp(n, 2*GF2E::WordLength());

      u = a;
      v = b;

      for (;;) {
         d0 = deg(u);
         d1 = deg(v);
         lc = LeadCoeff(v);

         PlainRem(u, u, v, tmp);
         swap(u, v);

         d2 = deg(v);
         if (d2 >= 0) {
            power(lc, lc, d0-d2);
            mul(res, res, lc);
            if (d0 & d1 & 1) negate(res, res);
         }
         else {
            if (d1 == 0) {
               power(lc, lc, d0);
               mul(res, res, lc);
            }
            else
               clear(res);
        
            break;
         }
      }

      rres = res;
   }
}

void resultant(GF2E& rres, const GF2EX& a, const GF2EX& b)
{
   PlainResultant(rres, a, b); 
}


void NormMod(GF2E& x, const GF2EX& a, const GF2EX& f)
{
   if (deg(f) <= 0 || deg(a) >= deg(f)) 
      Error("norm: bad args");

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

   GF2E t;
   resultant(t, f, a);
   if (!IsOne(LeadCoeff(f))) {
      GF2E t1;
      power(t1, LeadCoeff(f), deg(a));
      inv(t1, t1);
      mul(t, t, t1);
   }

   x = t;
}



// tower stuff...

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

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

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

      if (coeff(v, i) != 0) {
         for (j = 0; j < m; j++) {
            add(t[j], t[j], h[j]);
         }
      }
   }

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

   x.normalize();
}



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


   GF2EX s, t;
   vec_GF2E scratch;
   scratch.SetLength(deg(F));

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

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

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


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

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


   GF2EXArgument A;

   build(A, h, F, m);

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

void PrepareProjection(vec_vec_GF2& tt, const vec_GF2E& s,
                       const vec_GF2& proj)
{
   long l = s.length();
   tt.SetLength(l);

   GF2XTransMultiplier M;
   long i;

   for (i = 0; i < l; i++) {
      build(M, rep(s[i]), GF2E::modulus());
      UpdateMap(tt[i], proj, M, GF2E::modulus());
   }
}

void ProjectedInnerProduct(GF2& x, const vec_GF2E& a, 
                           const vec_vec_GF2& b)
{
   long n = min(a.length(), b.length());

   GF2 t, res;

   res = 0;

   long i;
   for (i = 0; i < n; i++) {
      project(t, b[i], rep(a[i]));
      res += t;
   }

   x = res;
}



void PrecomputeProj(vec_GF2& proj, const GF2X& f)
{
   long n = deg(f);

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

   if (ConstTerm(f) != 0) {
      proj.SetLength(1);
      proj[0] = 1;
   }
   else {
      proj.SetLength(n);
      clear(proj);
      proj[n-1] = 1;
   }
}

void ProjectPowersTower(vec_GF2& x, const vec_GF2E& a, long k,
                   const GF2EXArgument& H, const GF2EXModulus& F,
                   const vec_GF2& proj)

{
   long n = F.n;

   if (a.length() > n || k < 0) 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);

   vec_GF2E s(INIT_SIZE, n);
   s = a;

   x.SetLength(k);

   vec_vec_GF2 tt;

   for (long i = 0; i <= l; i++) {
      long m1 = min(m, k-i*m);

      PrepareProjection(tt, s, proj);

      for (long j = 0; j < m1; j++) {
         GF2 r;
         ProjectedInnerProduct(r, H.H[j].rep, tt);
         x.put(i*m + j, r);
      }
      if (i < l)
         UpdateMap(s, s, M, F);
   }
}




void ProjectPowersTower(vec_GF2& x, const vec_GF2E& a, long k,
                   const GF2EX& h, const GF2EXModulus& F,
                   const vec_GF2& proj)

{
   if (a.length() > F.n || k < 0) Error("ProjectPowers: bad args");

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

   long m = SqrRoot(k);

   GF2EXArgument H;

   build(H, h, F, m);
   ProjectPowersTower(x, a, k, H, F, proj);
}


void DoMinPolyTower(GF2X& h, const GF2EX& g, const GF2EXModulus& F, long m,
               const vec_GF2E& R, const vec_GF2& proj)
{
   vec_GF2 x;

   ProjectPowersTower(x, R, 2*m, g, F, proj);
   
   MinPolySeq(h, x, m);
}


void ProbMinPolyTower(GF2X& h, const GF2EX& g, const GF2EXModulus& F, 
                      long m)
{
   long n = F.n;
   if (m < 1 || m > n*GF2E::degree()) Error("ProbMinPoly: bad args");

   vec_GF2E R;
   R.SetLength(n);
   long i;
   for (i = 0; i < n; i++) random(R[i]);

   vec_GF2 proj;
   PrecomputeProj(proj, GF2E::modulus());

   DoMinPolyTower(h, g, F, m, R, proj);
}

void ProbMinPolyTower(GF2X& h, const GF2EX& g, const GF2EXModulus& F, 
                      long m, const vec_GF2& proj)
{
   long n = F.n;
   if (m < 1 || m > n*GF2E::degree()) Error("ProbMinPoly: bad args");

   vec_GF2E R;
   R.SetLength(n);
   long i;
   for (i = 0; i < n; i++) random(R[i]);

   DoMinPolyTower(h, g, F, m, R, proj);
}

void MinPolyTower(GF2X& hh, const GF2EX& g, const GF2EXModulus& F, long m)
{
   GF2X h;
   GF2EX h1;
   long n = F.n;
   if (m < 1 || m > n*GF2E::degree()) {
      Error("MinPoly: bad args");
   }

   vec_GF2 proj;
   PrecomputeProj(proj, GF2E::modulus());

   /* probabilistically compute min-poly */

   ProbMinPolyTower(h, g, F, m, proj);
   if (deg(h) == m) { hh = h; return; }
   CompTower(h1, h, g, F);
   if (IsZero(h1)) { hh = h; return; }

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

   long i;

   GF2X h2;
   GF2EX h3;
   vec_GF2E R;
   GF2EXTransMultiplier H1;
   

   for (;;) {
      R.SetLength(n);
      for (i = 0; i < n; i++) random(R[i]);
      build(H1, h1, F);
      UpdateMap(R, R, H1, F);
      DoMinPolyTower(h2, g, F, m-deg(h), R, proj);

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

void IrredPolyTower(GF2X& h, const GF2EX& g, const GF2EXModulus& F, long m)
{
   if (m < 1 || m > deg(F)*GF2E::degree()) Error("IrredPoly: bad args");

   vec_GF2E R;
   R.SetLength(1);
   R[0] = 1;

   vec_GF2 proj;
   proj.SetLength(1);
   proj.put(0, 1);

   DoMinPolyTower(h, g, F, m, R, proj);
}

NTL_END_IMPL