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 < 0) { clear(x); return; }
   if (NTL_OVERFLOW(hi, 1, 0)) 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