void NewtonInv(ZZ_pEX& c, const ZZ_pEX& a, long e)
{
   ZZ_pE x;

   inv(x, ConstTerm(a));

   if (e == 1) {
      conv(c, x);
      return;
   }

   static vec_long E;
   E.SetLength(0);
   append(E, e);
   while (e > 1) {
      e = (e+1)/2;
      append(E, e);
   }

   long L = E.length();

   ZZ_pEX g, g0, g1, g2;


   g.rep.SetMaxLength(E[0]);
   g0.rep.SetMaxLength(E[0]);
   g1.rep.SetMaxLength((3*E[0]+1)/2);
   g2.rep.SetMaxLength(E[0]);

   conv(g, x);

   long i;

   for (i = L-1; i > 0; i--) {
      // lift from E[i] to E[i-1]

      long k = E[i];
      long l = E[i-1]-E[i];

      trunc(g0, a, k+l);

      mul(g1, g0, g);
      RightShift(g1, g1, k);
      trunc(g1, g1, l);

      mul(g2, g1, g);
      trunc(g2, g2, l);
      LeftShift(g2, g2, k);

      sub(g, g, g2);
   }

   c = g;
}

void InvTrunc(ZZ_pEX& c, const ZZ_pEX& a, long e)
{
   if (e < 0) Error("InvTrunc: bad args");
   if (e == 0) {
      clear(c);
      return;
   }

   if (NTL_OVERFLOW(e, 1, 0))
      Error("overflow in InvTrunc");

   NewtonInv(c, a, e);
}




const long ZZ_pEX_MOD_PLAIN = 0;
const long ZZ_pEX_MOD_MUL = 1;


void build(ZZ_pEXModulus& F, const ZZ_pEX& f)
{
   long n = deg(f);

   if (n <= 0) Error("build(ZZ_pEXModulus,ZZ_pEX): deg(f) <= 0");

   if (NTL_OVERFLOW(n, ZZ_pE::degree(), 0))
      Error("build(ZZ_pEXModulus,ZZ_pEX): overflow");
   

   F.tracevec.SetLength(0);

   F.f = f;
   F.n = n;

   if (F.n < ZZ_pE::ModCross()) {
      F.method = ZZ_pEX_MOD_PLAIN;
   }
   else {
      F.method = ZZ_pEX_MOD_MUL;
      ZZ_pEX P1;
      ZZ_pEX P2;

      CopyReverse(P1, f, n);
      InvTrunc(P2, P1, n-1);
      CopyReverse(P1, P2, n-2);
      trunc(F.h0, P1, n-2);
      trunc(F.f0, f, n);
      F.hlc = ConstTerm(P2);
   }
}



ZZ_pEXModulus::ZZ_pEXModulus()
{
   n = -1;
   method = ZZ_pEX_MOD_PLAIN;
}


ZZ_pEXModulus::~ZZ_pEXModulus() 
{ 
}



ZZ_pEXModulus::ZZ_pEXModulus(const ZZ_pEX& ff)
{
   n = -1;
   method = ZZ_pEX_MOD_PLAIN;

   build(*this, ff);
}


void UseMulRem21(ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEXModulus& F)
{
   ZZ_pEX P1;
   ZZ_pEX P2;

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
   add(P2, P2, P1);
   mul(P1, P2, F.f0);
   trunc(P1, P1, F.n);
   trunc(r, a, F.n);
   sub(r, r, P1);
}

void UseMulDivRem21(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEXModulus& F)
{
   ZZ_pEX P1;
   ZZ_pEX P2;

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
   add(P2, P2, P1);
   mul(P1, P2, F.f0);
   trunc(P1, P1, F.n);
   trunc(r, a, F.n);
   sub(r, r, P1);
   q = P2;
}

void UseMulDiv21(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEXModulus& F)
{
   ZZ_pEX P1;
   ZZ_pEX P2;

   RightShift(P1, a, F.n);
   mul(P2, P1, F.h0);
   RightShift(P2, P2, F.n-2);
   if (!IsOne(F.hlc)) mul(P1, P1, F.hlc);
   add(P2, P2, P1);
   q = P2;

}


void rem(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEXModulus& F)
{
   if (F.method == ZZ_pEX_MOD_PLAIN) {
      PlainRem(x, a, F.f);
      return;
   }

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

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

   ZZ_pEX buf(INIT_SIZE, 2*n-1);

   long a_len = da+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();

      UseMulRem21(buf, buf, F);

      a_len -= amt;
   }

   x = buf;
}

void DivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEXModulus& F)
{
   if (F.method == ZZ_pEX_MOD_PLAIN) {
      PlainDivRem(q, r, a, F.f);
      return;
   }

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

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

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

   ZZ_pEX 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();

      UseMulDivRem21(qbuf, buf, buf, F);
      long dl = qbuf.rep.length();
      a_len = a_len - amt;
      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;
   }

   r = buf;

   qq.normalize();
   q = qq;
}

void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEXModulus& F)
{
   if (F.method == ZZ_pEX_MOD_PLAIN) {
      PlainDiv(q, a, F.f);
      return;
   }

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

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

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

   ZZ_pEX 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(ZZ_pEX& c, const ZZ_pEX& a, const ZZ_pEX& b, const ZZ_pEXModulus& F)
{
   if (deg(a) >= F.n || deg(b) >= F.n) Error("MulMod: bad args");

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


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

   ZZ_pEX t;
   sqr(t, a);
   rem(c, t, F);
}



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

void UseMulDivRem(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
{
   ZZ_pEX P1;
   ZZ_pEX 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);
   sub(P1, a, P1);
   
   r = P1;
   q = P2;
}

void UseMulDiv(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pEX& b)
{
   ZZ_pEX P1;
   ZZ_pEX 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(ZZ_pEX& q, ZZ_pEX& r, const ZZ_pEX& a, const ZZ_pEX& b)
{
   long sa = a.rep.length();
   long sb = b.rep.length();

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

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

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

void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_pE& b)
{
   ZZ_pE T;
   inv(T, b);
   mul(q, a, T);
}

void div(ZZ_pEX& q, const ZZ_pEX& a, const ZZ_p& b)
{
   NTL_ZZ_pRegister(T);
   inv(T, b);
   mul(q, a, T);
}

void div(ZZ_pEX& q, const ZZ_pEX& a, long b)
{
   NTL_ZZ_pRegister(T);
   T = b;
   inv(T, T);
   mul(q, a, T);
}

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

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

void GCD(ZZ_pEX& x, const ZZ_pEX& a, const ZZ_pEX& b)
{
   ZZ_pE t;

   if (IsZero(b))
      x = a;
   else if (IsZero(a))
      x = b;
   else {
      long n = max(deg(a),deg(b)) + 1;
      ZZ_pEX u(INIT_SIZE, n), v(INIT_SIZE, n);

      vec_ZZ_pX tmp;
      SetSize(tmp, n, 2*ZZ_pE::degree());

      u = a;
      v = b;
      do {
         PlainRem(u, u, v, tmp);
         swap(u, v);
      } while (!IsZero(v));

      x = u;
   }

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

   /* make gcd monic */


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



         

void XGCD(ZZ_pEX& d, ZZ_pEX& s, ZZ_pEX& t, const ZZ_pEX& a, const ZZ_pEX& b)
{
   ZZ_pE z;


   if (IsZero(b)) {
      set(s);
      clear(t);
      d = a;
   }
   else if (IsZero(a)) {
      clear(s);
      set(t);
      d = b;
   }
   else {
      long e = max(deg(a), deg(b)) + 1;

      ZZ_pEX temp(INIT_SIZE, e), u(INIT_SIZE, e), v(INIT_SIZE, e), 
            u0(INIT_SIZE, e), v0(INIT_SIZE, e), 
            u1(INIT_SIZE, e), v1(INIT_SIZE, e), 
            u2(INIT_SIZE, e), v2(INIT_SIZE, e), q(INIT_SIZE, e);


      set(u1); clear(v1);
      clear(u2); set(v2);
      u = a; v = b;

      do {
         DivRem(q, u, u, v);
         swap(u, v);
         u0 = u2;
         v0 = v2;
         mul(temp, q, u2);
         sub(u2, u1, temp);
         mul(temp, q, v2);
         sub(v2, v1, temp);
         u1 = u0;
         v1 = v0;
      } while (!IsZero(v));

      d = u;
      s = u1;
      t = v1;
   }

   if (IsZero(d)) return;
   if (IsOne(LeadCoeff(d))) return;

   /* make gcd monic */

   inv(z, LeadCoeff(d));
   mul(d, d, z);
   mul(s, s, z);
   mul(t, t, z);
}

NTL_vector_impl(ZZ_pEX,vec_ZZ_pEX)

NTL_eq_vector_impl(ZZ_pEX,vec_ZZ_pEX)

NTL_io_vector_impl(ZZ_pEX,vec_ZZ_pEX)

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

   if (n <= 0) return;

   negate(a[0], a[0]);

   for (k = 1; k <= n-1; k++) {
      negate(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(ZZ_pEX& x, const vec_ZZ_pE& 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(ZZ_pE& b, const ZZ_pEX& f, const ZZ_pE& a)
// does a Horner evaluation
{
   ZZ_pE 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_ZZ_pE& b, const ZZ_pEX& f, const vec_ZZ_pE& a)
// naive algorithm:  repeats Horner
{
   if (&b == &f.rep) {
      vec_ZZ_pE 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(ZZ_pEX& f, const vec_ZZ_pE& a, const vec_ZZ_pE& b)
{
   long m = a.length();
   if (b.length() != m) Error("interpolate: vector length mismatch");

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

   vec_ZZ_pE prod;
   prod = a;

   ZZ_pE t1, t2;