#include <NTL/lzz_pX.h>

#include <NTL/new.h>

NTL_START_IMPL



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

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

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



void zz_pXMatrix::operator=(const zz_pXMatrix& M)
{
   elts[0][0] = M.elts[0][0];
   elts[0][1] = M.elts[0][1];
   elts[1][0] = M.elts[1][0];
   elts[1][1] = M.elts[1][1];
}


void RightShift(zz_pX& x, const zz_pX& a, long n)
{
   if (IsZero(a)) {
      clear(x);
      return;
   }

   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(zz_pX& x, const zz_pX& a, long n)
{
   if (IsZero(a)) {
      clear(x);
      return;
   }

   if (n < 0) {
      if (n < -NTL_MAX_LONG) 
         clear(x);
      else
         RightShift(x, a, -n);
      return;
   }

   if (NTL_OVERFLOW(n, 1, 0))
      Error("overflow in LeftShift");

   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(zz_pX& U, const zz_pX& 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();
}

void ShiftSub(zz_pX& U, const zz_pX& 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++)
      sub(U.rep[i+n], U.rep[i+n], V.rep[i]);

   U.normalize();
}

void mul(zz_pX& U, zz_pX& V, const zz_pXMatrix& M)
// (U, V)^T = M*(U, V)^T
{
   long d = deg(U) - deg(M(1,1));
   long k = NextPowerOfTwo(d - 1);

   // When the GCD algorithm is run on polynomials of degree n, n-1, 
   // where n is a power of two, then d-1 is likely to be a power of two.
   // It would be more natural to set k = NextPowerOfTwo(d+1), but this
   // would be much less efficient in this case.

   long n = (1L << k);
   long xx;
   zz_p a0, a1, b0, b1, c0, d0, u0, u1, v0, v1, nu0, nu1, nv0;
   zz_p t1, t2;

   if (n == d-1)
      xx = 1;
   else if (n == d)
      xx = 2;
   else 
      xx = 3;

   switch (xx) {
   case 1:
      GetCoeff(a0, M(0,0), 0);
      GetCoeff(a1, M(0,0), 1);
      GetCoeff(b0, M(0,1), 0);
      GetCoeff(b1, M(0,1), 1);
      GetCoeff(c0, M(1,0), 0);
      GetCoeff(d0, M(1,1), 0);

      GetCoeff(u0, U, 0);
      GetCoeff(u1, U, 1);
      GetCoeff(v0, V, 0);
      GetCoeff(v1, V, 1);

      mul(t1, (a0), (u0));
      mul(t2, (b0), (v0));
      add(t1, t1, t2); 
      nu0 = t1;

      mul(t1, (a1), (u0));
      mul(t2, (a0), (u1));
      add(t1, t1, t2);
      mul(t2, (b1), (v0));
      add(t1, t1, t2);
      mul(t2, (b0), (v1));
      add(t1, t1, t2);
      nu1 = t1;

      mul(t1, (c0), (u0));
      mul(t2, (d0), (v0));
      add (t1, t1, t2);
      nv0 = t1;
   
      break;

   case 2:
      GetCoeff(a0, M(0,0), 0);
      GetCoeff(b0, M(0,1), 0);

      GetCoeff(u0, U, 0);
      GetCoeff(v0, V, 0);

      mul(t1, (a0), (u0));
      mul(t2, (b0), (v0));
      add(t1, t1, t2); 
      nu0 = t1;

      break;

   case 3:
      break;

   }

   fftRep RU(INIT_SIZE, k), RV(INIT_SIZE, k), R1(INIT_SIZE, k), 
          R2(INIT_SIZE, k);

   TofftRep(RU, U, k);  
   TofftRep(RV, V, k);  

   TofftRep(R1, M(0,0), k);
   mul(R1, R1, RU);
   TofftRep(R2, M(0,1), k);
   mul(R2, R2, RV);
   add(R1, R1, R2);
   FromfftRep(U, R1, 0, d);

   TofftRep(R1, M(1,0), k);
   mul(R1, R1, RU);
   TofftRep(R2, M(1,1), k);
   mul(R2, R2, RV);
   add(R1, R1, R2);
   FromfftRep(V, R1, 0, d-1);

   // now fix-up results

   switch (xx) {
   case 1:
      GetCoeff(u0, U, 0);
      sub(u0, u0, nu0);
      SetCoeff(U, d-1, u0);
      SetCoeff(U, 0, nu0);

      GetCoeff(u1, U, 1);
      sub(u1, u1, nu1);
      SetCoeff(U, d, u1);
      SetCoeff(U, 1, nu1);

      GetCoeff(v0, V, 0);
      sub(v0, v0, nv0);
      SetCoeff(V, d-1, v0);
      SetCoeff(V, 0, nv0);

      break;
      

   case 2:
      GetCoeff(u0, U, 0);
      sub(u0, u0, nu0);
      SetCoeff(U, d, u0);
      SetCoeff(U, 0, nu0);

      break;

   }
}


void mul(zz_pXMatrix& A, zz_pXMatrix& B, zz_pXMatrix& C)
// A = B*C, B and C are destroyed
{
   long db = deg(B(1,1));
   long dc = deg(C(1,1));
   long da = db + dc;

   long k = NextPowerOfTwo(da+1);

   fftRep B00, B01, B10, B11, C0, C1, T1, T2;
   
   TofftRep(B00, B(0,0), k); B(0,0).kill();
   TofftRep(B01, B(0,1), k); B(0,1).kill();
   TofftRep(B10, B(1,0), k); B(1,0).kill();
   TofftRep(B11, B(1,1), k); B(1,1).kill();

   TofftRep(C0, C(0,0), k);  C(0,0).kill();
   TofftRep(C1, C(1,0), k);  C(1,0).kill();

   mul(T1, B00, C0);
   mul(T2, B01, C1);
   add(T1, T1, T2);
   FromfftRep(A(0,0), T1, 0, da);

   mul(T1, B10, C0);
   mul(T2, B11, C1);
   add(T1, T1, T2);
   FromfftRep(A(1,0), T1, 0, da);

   TofftRep(C0, C(0,1), k);  C(0,1).kill();
   TofftRep(C1, C(1,1), k);  C(1,1).kill();

   mul(T1, B00, C0);
   mul(T2, B01, C1);
   add(T1, T1, T2);
   FromfftRep(A(0,1), T1, 0, da);

   mul(T1, B10, C0);
   mul(T2, B11, C1);
   add(T1, T1, T2);
   FromfftRep(A(1,1), T1, 0, da);
}

void IterHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red)
{
   M_out(0,0).SetMaxLength(d_red);
   M_out(0,1).SetMaxLength(d_red);
   M_out(1,0).SetMaxLength(d_red);
   M_out(1,1).SetMaxLength(d_red);

   set(M_out(0,0));   clear(M_out(0,1));
   clear(M_out(1,0)); set(M_out(1,1));

   long goal = deg(U) - d_red;

   if (deg(V) <= goal)
      return;

   zz_pX Q, t(INIT_SIZE, d_red);

   while (deg(V) > goal) {
      PlainDivRem(Q, U, U, V);
      swap(U, V);

      mul(t, Q, M_out(1,0));
      sub(t, M_out(0,0), t);
      M_out(0,0) = M_out(1,0);
      M_out(1,0) = t;

      mul(t, Q, M_out(1,1));
      sub(t, M_out(0,1), t);
      M_out(0,1) = M_out(1,1);
      M_out(1,1) = t;
   }
}
   


void HalfGCD(zz_pXMatrix& M_out, const zz_pX& U, const zz_pX& V, long d_red)
{
   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      set(M_out(0,0));   clear(M_out(0,1));
      clear(M_out(1,0)); set(M_out(1,1));
 
      return;
   }


   long n = deg(U) - 2*d_red + 2;
   if (n < 0) n = 0;

   zz_pX U1, V1;

   RightShift(U1, U, n);
   RightShift(V1, V, n);

   if (d_red <= NTL_zz_pX_HalfGCD_CROSSOVER) {
      IterHalfGCD(M_out, U1, V1, d_red);
      return;
   }

   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   zz_pXMatrix M1;

   HalfGCD(M1, U1, V1, d1);
   mul(U1, V1, M1);

   long d2 = deg(V1) - deg(U) + n + d_red;

   if (IsZero(V1) || d2 <= 0) {
      M_out = M1;
      return;
   }


   zz_pX Q;
   zz_pXMatrix M2;

   DivRem(Q, U1, U1, V1);
   swap(U1, V1);

   HalfGCD(M2, U1, V1, d2);

   zz_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,0));
   sub(t, M1(0,0), t);
   swap(M1(0,0), M1(1,0));
   swap(M1(1,0), t);

   t.kill();

   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,1));
   sub(t, M1(0,1), t);
   swap(M1(0,1), M1(1,1));
   swap(M1(1,1), t);

   t.kill();

   mul(M_out, M2, M1); 
}




void XHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red)
{
   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      set(M_out(0,0));   clear(M_out(0,1));
      clear(M_out(1,0)); set(M_out(1,1));
 
      return;
   }

   long du = deg(U);

   if (d_red <= NTL_zz_pX_HalfGCD_CROSSOVER) {
      IterHalfGCD(M_out, U, V, d_red);
      return;
   }

   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   zz_pXMatrix M1;

   HalfGCD(M1, U, V, d1);
   mul(U, V, M1);

   long d2 = deg(V) - du + d_red;

   if (IsZero(V) || d2 <= 0) {
      M_out = M1;
      return;
   }


   zz_pX Q;
   zz_pXMatrix M2;

   DivRem(Q, U, U, V);
   swap(U, V);

   XHalfGCD(M2, U, V, d2);

   zz_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,0));
   sub(t, M1(0,0), t);
   swap(M1(0,0), M1(1,0));
   swap(M1(1,0), t);

   t.kill();

   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,1));
   sub(t, M1(0,1), t);
   swap(M1(0,1), M1(1,1));
   swap(M1(1,1), t);

   t.kill();

   mul(M_out, M2, M1); 
}

void HalfGCD(zz_pX& U, zz_pX& V)
{
   long d_red = (deg(U)+1)/2;

   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      return;
   }

   long du = deg(U);


   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   zz_pXMatrix M1;

   HalfGCD(M1, U, V, d1);
   mul(U, V, M1);

   long d2 = deg(V) - du + d_red;

   if (IsZero(V) || d2 <= 0) {
      return;
   }

   M1(0,0).kill();
   M1(0,1).kill();
   M1(1,0).kill();
   M1(1,1).kill();


   zz_pX Q;

   DivRem(Q, U, U, V);
   swap(U, V);

   HalfGCD(M1, U, V, d2);

   mul(U, V, M1); 
}


void GCD(zz_pX& d, const zz_pX& u, const zz_pX& v)
{
   zz_pX u1, v1;

   u1 = u;
   v1 = v;

   if (deg(u1) == deg(v1)) {
      if (IsZero(u1)) {
         clear(d);
         return;
      }

      rem(v1, v1, u1);
   }
   else if (deg(u1) < deg(v1)) {
      swap(u1, v1);
   }

   // deg(u1) > deg(v1)

   while (deg(u1) > NTL_zz_pX_GCD_CROSSOVER && !IsZero(v1)) {
      HalfGCD(u1, v1);

      if (!IsZero(v1)) {
         rem(u1, u1, v1);
         swap(u1, v1);
      }
   }

   PlainGCD(d, u1, v1);
}



void XGCD(zz_pX& d, zz_pX& s, zz_pX& t, const zz_pX& a, const zz_pX& b)
{
   zz_p w;

   if (IsZero(a) && IsZero(b)) {
      clear(d);
      set(s);
      clear(t);
      return;
   }

   zz_pX U, V, Q;

   U = a;
   V = b;

   long flag = 0;

   if (deg(U) == deg(V)) {
      DivRem(Q, U, U, V);
      swap(U, V);
      flag = 1;
   }
   else if (deg(U) < deg(V)) {
      swap(U, V);
      flag = 2;
   }

   zz_pXMatrix M;

   XHalfGCD(M, U, V, deg(U)+1);

   d = U;

   if (flag == 0) {
      s = M(0,0); 
      t = M(0,1);
   }
   else if (flag == 1) {
      s = M(0,1);
      mul(t, Q, M(0,1));
      sub(t, M(0,0), t);
   }
   else {  /* flag == 2 */
      s = M(0,1);
      t = M(0,0);
   }

   // normalize

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

      





void IterBuild(zz_p* a, long n)
{
   long i, k;
   zz_p 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 mul(zz_p* x, const zz_p* a, const zz_p* b, long n)
{
   zz_p t, accum;

   long i, j, jmin, jmax;

   long d = 2*n-1;

   for (i = 0; i <= d; i++) {
      jmin = max(0, i-(n-1));
      jmax = min(n-1, i);
      clear(accum);
      for (j = jmin; j <= jmax; j++) {
         mul(t, (a[j]), (b[i-j]));
         add(accum, accum, t);
      }
      if (i >= n) {
         add(accum, accum, (a[i-n]));
         add(accum, accum, (b[i-n]));
      }

      x[i] = accum;
   }
}


void BuildFromRoots(zz_pX& x, const vec_zz_p& a)
{
   long n = a.length();

   if (n == 0) {
      set(x);
      return;
   }

   long k0 = NextPowerOfTwo(NTL_zz_pX_MUL_CROSSOVER)-1;
   long crossover = 1L << k0;

   if (n <= NTL_zz_pX_MUL_CROSSOVER) {
      x.rep.SetMaxLength(n+1);
      x.rep = a;
      IterBuild(&x.rep[0], n);
      x.rep.SetLength(n+1);
      SetCoeff(x, n);
      return;
   }

   long k = NextPowerOfTwo(n);

   long m = 1L << k;
   long i, j;
   long l, width;

   zz_pX b(INIT_SIZE, m+1);

   b.rep = a;
   b.rep.SetLength(m+1);
   for (i = n; i < m; i++)
      clear(b.rep[i]);

   set(b.rep[m]);
   
   fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);


   zz_p t1, one;
   set(one);

   vec_zz_p G(INIT_SIZE, crossover), H(INIT_SIZE, crossover);
   zz_p *g = G.elts();
   zz_p *h = H.elts();
   zz_p *tmp;
   
   for (i = 0; i < m; i+= crossover) {
      for (j = 0; j < crossover; j++)
         negate(g[j], b.rep[i+j]);

      if (k0 > 0) {
         for (j = 0; j < crossover; j+=2) {
            mul(t1, g[j], g[j+1]);
            add(g[j+1], g[j], g[j+1]);
            g[j] = t1;
         }
      }
   
      for (l = 1; l < k0; l++) {
         width = 1L << l;

         for (j = 0; j < crossover; j += 2*width)
            mul(&h[j], &g[j], &g[j+width], width);
      
         tmp = g; g = h; h = tmp;
      }

      for (j = 0; j < crossover; j++)
         b.rep[i+j] = g[j];
   }

   for (l = k0; l < k; l++) {
      width = 1L << l;
      for (i = 0; i < m; i += 2*width) {
         t1 = b.rep[i+width];
         set(b.rep[i+width]);
         TofftRep(R1, b, l+1, i, i+width);
         b.rep[i+width] = t1;
         t1 = b.rep[i+2*width];
         set(b.rep[i+2*width]);
         TofftRep(R2, b, l+1, i+width, i+2*width);
         b.rep[i+2*width] = t1;
         mul(R1, R1, R2);
         FromfftRep(&b.rep[i], R1, 0, 2*width-1);
         sub(b.rep[i], b.rep[i], one);
      }
   }

   x.rep.SetLength(n+1);
   long delta = m-n;
   for (i = 0; i <= n; i++)
     x.rep[i] = b.rep[i+delta];

   // no need to normalize
}



void eval(zz_p& b, const zz_pX& f, zz_p a)
// does a Horner evaluation
{
   zz_p 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_p& b, const zz_pX& f, const vec_zz_p& a)
// naive algorithm:  repeats Horner
{
   if (&b == &f.rep) {
      vec_zz_p 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_pX& f, const vec_zz_p& a, const vec_zz_p& b)
{
   long m = a.length();
   if (b.length() != m) Error("interpolate: vector length mismatch");

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

   vec_zz_p prod;
   prod = a;

   zz_p t1, t2;

   long k, i;

   vec_zz_p res;
   res.SetLength(m);

   for (k = 0; k < m; k++) {

      const zz_p& 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;
}

NTL_vector_impl(zz_pX,vec_zz_pX)

NTL_eq_vector_impl(zz_pX,vec_zz_pX)



   
void InnerProduct(zz_pX& x, const vec_zz_p& v, long low, long high, 
                   const vec_zz_pX& H, long n, vec_zz_p& t)
{
   zz_p s;
   long i, j;

   zz_p *tp = t.elts();

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


   long p = zz_p::modulus();
   double pinv = zz_p::ModulusInverse();

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

      long W = rep(w);
      mulmod_precon_t Wpinv = PrepMulModPrecon(W, p, pinv); // ((double) W)*pinv;
      const zz_p *hp = h.elts();

      for (j = 0; j < m; j++) {
         long S = MulModPrecon(rep(hp[j]), W, p, Wpinv);
         S = AddMod(S, rep(tp[j]), p);
         tp[j].LoopHole() = S;
      }
   }

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


void CompMod(zz_pX& x, const zz_pX& g, const zz_pXArgument& A, 
             const zz_pXModulus& F)
{
   if (deg(g) <= 0) {
      x = g;
      return;
   }


   zz_pX s, t;
   vec_zz_p scratch(INIT_SIZE, F.n);

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

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

   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(zz_pXArgument& A, const zz_pX& h, const zz_pXModulus& F, long m)
{
   if (m <= 0 || deg(h) >= F.n) Error("build: bad args");

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

   long i;

   if (zz_pXArgBound > 0) {
      double sz = 1;
      sz = sz*F.n;
      sz = sz+6;
      sz = sz*(sizeof (long));
      sz = sz/1024;
      m = min(m, long(zz_pXArgBound/sz));
      m = max(m, 1);
   }

   zz_pXMultiplier M;

   build(M, h, F);

   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], M, F);
}




long zz_pXArgBound = 0;


void CompMod(zz_pX& x, const zz_pX& g, const zz_pX& h, const zz_pXModulus& F)
   // x = g(h) mod f
{
   long m = SqrRoot(g.rep.length());

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

   zz_pXArgument A;

   build(A, h, F, m);

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




void Comp2Mod(zz_pX& x1, zz_pX& x2, const zz_pX& g1, const zz_pX& g2,
              const zz_pX& h, const zz_pXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length());

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

   zz_pXArgument A;

   build(A, h, F, m);

   zz_pX xx1, xx2;

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

   x1 = xx1;
   x2 = xx2;
}

void Comp3Mod(zz_pX& x1, zz_pX& x2, zz_pX& x3, 
              const zz_pX& g1, const zz_pX& g2, const zz_pX& g3,
              const zz_pX& h, const zz_pXModulus& F)

{
   long m = SqrRoot(g1.rep.length() + g2.rep.length() + g3.rep.length());

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

   zz_pXArgument A;

   build(A, h, F, m);

   zz_pX xx1, xx2, xx3;

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

   x1 = xx1;
   x2 = xx2;
   x3 = xx3;
}

static void StripZeroes(vec_zz_p& x)
{
   long n = x.length();
   while (n > 0 && IsZero(x[n-1]))
      n--;
   x.SetLength(n);
}


void PlainUpdateMap(vec_zz_p& xx, const vec_zz_p& a, 
                    const zz_pX& b, const zz_pX& f)
{
   long n = deg(f);
   long i, m;

   if (IsZero(b)) {
      xx.SetLength(0);
      return;
   }

   m = n-1 - deg(b);

   vec_zz_p x(INIT_SIZE, n);

   for (i = 0; i <= m; i++)
      InnerProduct(x[i], a, b.rep, i);

   if (deg(b) != 0) {
      zz_pX c(INIT_SIZE, n);
      LeftShift(c, b, m);

      for (i = m+1; i < n; i++) {
         MulByXMod(c, c, f);
         InnerProduct(x[i], a, c.rep);
      }
   }

   xx = x;
}
   



void UpdateMap(vec_zz_p& x, const vec_zz_p& aa, 
               const zz_pXMultiplier& B, const zz_pXModulus& F)
{
   long n = F.n;

   vec_zz_p a;
   a = aa;
   StripZeroes(a);

   if (a.length() > n) Error("UpdateMap: bad args");
   long i;

   if (!B.UseFFT) {
      PlainUpdateMap(x, a, B.b, F.f);
      StripZeroes(x);
      return;
   }

   fftRep R1(INIT_SIZE, F.k), R2(INIT_SIZE, F.l);
   vec_zz_p V1(INIT_SIZE, n);


   RevTofftRep(R1, a, F.k, 0, a.length()-1, 0);
   mul(R2, R1, F.FRep);
   RevFromfftRep(V1, R2, 0, n-2);
   for (i = 0; i <= n-2; i++)  negate(V1[i], V1[i]);
   RevTofftRep(R2, V1, F.l, 0, n-2, n-1);
   mul(R2, R2, B.B1);
   mul(R1, R1, B.B2);

   AddExpand(R2, R1);
   RevFromfftRep(x, R2, 0, n-1);
   StripZeroes(x);
}

   

void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
                   const zz_pXArgument& H, const zz_pXModulus& F)

{
   long n = F.n;

   if (a.length() > n || k < 0 || NTL_OVERFLOW(k, 1, 0))
      Error("ProjectPowers: bad args");

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

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

   vec_zz_p s(INIT_SIZE, n);
   s = a;
   StripZeroes(s);

   x.SetLength(k);

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



void ProjectPowers(vec_zz_p& x, const vec_zz_p& a, long k,
                   const zz_pX& h, const zz_pXModulus& F)

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

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

   long m = SqrRoot(k);

   zz_pXArgument H;

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


void BerlekampMassey(zz_pX& h, const vec_zz_p& a, long m)
{
   zz_pX Lambda, Sigma, Temp;
   long L;
   zz_p Delta, Delta1, t1;
   long shamt;

   // 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(Delta1);
      dl = deg(Lambda);
      for (i = 0; i <= dl; i++) {
         mul(t1, Lambda.rep[i], a[r-i-1]);
         add(Delta1, Delta1, t1);
      }

      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;
         ShiftSub(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);
         ShiftSub(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 GCDMinPolySeq(zz_pX& h, const vec_zz_p& x, long m)
{
   long i;
   zz_pX a, b;
   zz_pXMatrix M;
   zz_p t;

   a.rep.SetLength(2*m);
   for (i = 0; i < 2*m; i++) a.rep[i] = x[2*m-1-i];
   a.normalize();

   SetCoeff(b, 2*m);

   HalfGCD(M, b, a, m+1);

   /* make monic */

   inv(t, LeadCoeff(M(1,1)));
   mul(h, M(1,1), t);
}


void MinPolySeq(zz_pX& h, const vec_zz_p& a, long m)
{
   if (m < 0 || NTL_OVERFLOW(m, 1, 0)) Error("MinPoly: bad args");
   if (a.length() < 2*m) Error("MinPoly: sequence too short");

   if (m > NTL_zz_pX_BERMASS_CROSSOVER)
      GCDMinPolySeq(h, a, m);
   else
      BerlekampMassey(h, a, m);
}


void DoMinPolyMod(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m,
               const vec_zz_p& R) 
{
   vec_zz_p x;

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


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

   long i;
   vec_zz_p R(INIT_SIZE, n);

   for (i = 0; i < n; i++) random(R[i]);
   DoMinPolyMod(h, g, F, m, R);
}

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

   long i;

   zz_pX h2, h3;
   zz_pXMultiplier H1;
   vec_zz_p R(INIT_SIZE, n);

   for (;;) {
      R.SetLength(n);
      for (i = 0; i < n; i++) random(R[i]);
      build(H1, h1, F);
      UpdateMap(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(zz_pX& h, const zz_pX& g, const zz_pXModulus& F, long m)
{
   vec_zz_p R(INIT_SIZE, 1);
   if (m < 1 || m > F.n) Error("IrredPoly: bad args");

   set(R[0]);
   DoMinPolyMod(h, g, F, m, R);
}



void diff(zz_pX& x, const zz_pX& 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++) {
      mul(x.rep[i], a.rep[i+1], i+1);
   }

   if (&x == &a)
      x.rep.SetLength(n);

   x.normalize();
}

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

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

   zz_p t;

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




      
void PlainMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n)
{
   zz_pX y;
   mul(y, a, b);
   trunc(x, y, n);
}


void FFTMulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n)
{
   if (IsZero(a) || IsZero(b)) {
      clear(x);
      return;
   }

   long d = deg(a) + deg(b);
   if (n > d + 1)
      n = d + 1;

   long k = NextPowerOfTwo(d + 1);
   fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);

   TofftRep(R1, a, k);
   TofftRep(R2, b, k);
   mul(R1, R1, R2);
   FromfftRep(x, R1, 0, n-1);
}

void MulTrunc(zz_pX& x, const zz_pX& a, const zz_pX& b, long n)
{
   if (n < 0) Error("MulTrunc: bad args");

   if (deg(a) <= NTL_zz_pX_MUL_CROSSOVER || deg(b) <= NTL_zz_pX_MUL_CROSSOVER)
      PlainMulTrunc(x, a, b, n);
   else
      FFTMulTrunc(x, a, b, n);
}

void PlainSqrTrunc(zz_pX& x, const zz_pX& a, long n)
{
   zz_pX y;
   sqr(y, a);
   trunc(x, y, n);
}


void FFTSqrTrunc(zz_pX& x, const zz_pX& a, long n)
{
   if (IsZero(a)) {
      clear(x);
      return;
   }

   long d = 2*deg(a);
   if (n > d + 1)
      n = d + 1;

   long k = NextPowerOfTwo(d + 1);
   fftRep R1(INIT_SIZE, k);

   TofftRep(R1, a, k);
   mul(R1, R1, R1);
   FromfftRep(x, R1, 0, n-1);
}

void SqrTrunc(zz_pX& x, const zz_pX& a, long n)
{
   if (n < 0) Error("SqrTrunc: bad args");

   if (deg(a) <= NTL_zz_pX_MUL_CROSSOVER)
      PlainSqrTrunc(x, a, n);
   else
      FFTSqrTrunc(x, a, n);
}



void FastTraceVec(vec_zz_p& S, const zz_pX& f)
{
   long n = deg(f);

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

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

   if (n == 1) {
      S.SetLength(1);
      set(S[0]);
      return;
   }
   
   long i;
   zz_pX f1;

   f1.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      f1.rep[i] = f.rep[n-i];
   f1.normalize();

   zz_pX f2;
   f2.rep.SetLength(n-1);
   for (i = 0; i <= n-2; i++)
      mul(f2.rep[i], f.rep[n-1-i], i+1);
   f2.normalize();

   zz_pX f3;
   InvTrunc(f3, f1, n-1);
   MulTrunc(f3, f3, f2, n-1);

   S.SetLength(n);

   S[0] = n;
   for (i = 1; i < n; i++)
      negate(S[i], coeff(f3, i-1));
}


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

   zz_pX f;
   f = ff;

   MakeMonic(f);

   long n = deg(f);

   S.SetLength(n);

   if (n == 0)
      return;

   long k, i;
   zz_p acc, t;

   const zz_p *fp = f.rep.elts();;
   zz_p *sp = S.elts();

   sp[0] = n;

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

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

      negate(sp[k], acc);
   }
}

void TraceVec(vec_zz_p& S, const zz_pX& f)
{
   if (deg(f) <= NTL_zz_pX_TRACE_CROSSOVER)
      PlainTraceVec(S, f);
   else
      FastTraceVec(S, f);
}

void ComputeTraceVec(const zz_pXModulus& F)
{
   vec_zz_p& S = *((vec_zz_p *) &F.tracevec);

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

   if (!F.UseFFT) {
      PlainTraceVec(S, F.f);
      return;
   }

   long i;
   long n = F.n;

   fftRep R;
   zz_pX P, g;

   g.rep.SetLength(n-1);
   for (i = 1; i < n; i++)
      mul(g.rep[n-i-1], F.f.rep[n-i], i); 
   g.normalize();

   TofftRep(R, g, F.l);
   mul(R, R, F.HRep);
   FromfftRep(P, R, n-2, 2*n-4);

   S.SetLength(n);

   S[0] = n;
   for (i = 1; i < n; i++)
      negate(S[i], coeff(P, n-1-i));
}

void TraceMod(zz_p& x, const zz_pX& a, const zz_pXModulus& F)
{
   long n = F.n;

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

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

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


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

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


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

      long n = max(deg(a),deg(b)) + 1;
      zz_pX u(INIT_SIZE, n), v(INIT_SIZE, n);

      u = a;
      v = b;

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

         PlainRem(u, u, v);
         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 ResIterHalfGCD(zz_pXMatrix& M_out, zz_pX& U, zz_pX& V, long d_red,
                    vec_zz_p& cvec, vec_long& dvec)
{
   M_out(0,0).SetMaxLength(d_red);
   M_out(0,1).SetMaxLength(d_red);
   M_out(1,0).SetMaxLength(d_red);
   M_out(1,1).SetMaxLength(d_red);

   set(M_out(0,0));   clear(M_out(0,1));
   clear(M_out(1,0)); set(M_out(1,1));

   long goal = deg(U) - d_red;

   if (deg(V) <= goal)
      return;

   zz_pX Q, t(INIT_SIZE, d_red);


   while (deg(V) > goal) {
      append(cvec, LeadCoeff(V));
      append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
      PlainDivRem(Q, U, U, V);
      swap(U, V);

      mul(t, Q, M_out(1,0));
      sub(t, M_out(0,0), t);
      M_out(0,0) = M_out(1,0);
      M_out(1,0) = t;

      mul(t, Q, M_out(1,1));
      sub(t, M_out(0,1), t);
      M_out(0,1) = M_out(1,1);
      M_out(1,1) = t;
   }
}
   


void ResHalfGCD(zz_pXMatrix& M_out, const zz_pX& U, const zz_pX& V, long d_red,
                vec_zz_p& cvec, vec_long& dvec)
{
   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      set(M_out(0,0));   clear(M_out(0,1));
      clear(M_out(1,0)); set(M_out(1,1));
 
      return;
   }


   long n = deg(U) - 2*d_red + 2;
   if (n < 0) n = 0;

   zz_pX U1, V1;

   RightShift(U1, U, n);
   RightShift(V1, V, n);

   if (d_red <= NTL_zz_pX_HalfGCD_CROSSOVER) { 
      ResIterHalfGCD(M_out, U1, V1, d_red, cvec, dvec);
      return;
   }

   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   zz_pXMatrix M1;

   ResHalfGCD(M1, U1, V1, d1, cvec, dvec);
   mul(U1, V1, M1);

   long d2 = deg(V1) - deg(U) + n + d_red;

   if (IsZero(V1) || d2 <= 0) {
      M_out = M1;
      return;
   }


   zz_pX Q;
   zz_pXMatrix M2;

   append(cvec, LeadCoeff(V1));
   append(dvec, dvec[dvec.length()-1]-deg(U1)+deg(V1));
   DivRem(Q, U1, U1, V1);
   swap(U1, V1);

   ResHalfGCD(M2, U1, V1, d2, cvec, dvec);

   zz_pX t(INIT_SIZE, deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,0));
   sub(t, M1(0,0), t);
   swap(M1(0,0), M1(1,0));
   swap(M1(1,0), t);

   t.kill();

   t.SetMaxLength(deg(M1(1,1))+deg(Q)+1);

   mul(t, Q, M1(1,1));
   sub(t, M1(0,1), t);
   swap(M1(0,1), M1(1,1));
   swap(M1(1,1), t);

   t.kill();

   mul(M_out, M2, M1); 
}

void ResHalfGCD(zz_pX& U, zz_pX& V, vec_zz_p& cvec, vec_long& dvec)
{
   long d_red = (deg(U)+1)/2;

   if (IsZero(V) || deg(V) <= deg(U) - d_red) {
      return;
   }

   long du = deg(U);


   long d1 = (d_red + 1)/2;
   if (d1 < 1) d1 = 1;
   if (d1 >= d_red) d1 = d_red - 1;

   zz_pXMatrix M1;

   ResHalfGCD(M1, U, V, d1, cvec, dvec);
   mul(U, V, M1);

   long d2 = deg(V) - du + d_red;

   if (IsZero(V) || d2 <= 0) {
      return;
   }

   M1(0,0).kill();
   M1(0,1).kill();
   M1(1,0).kill();
   M1(1,1).kill();


   zz_pX Q;

   append(cvec, LeadCoeff(V));
   append(dvec, dvec[dvec.length()-1]-deg(U)+deg(V));
   DivRem(Q, U, U, V);
   swap(U, V);

   ResHalfGCD(M1, U, V, d2, cvec, dvec);

   mul(U, V, M1); 
}


void resultant(zz_p& rres, const zz_pX& u, const zz_pX& v)
{
   if (deg(u) <= NTL_zz_pX_GCD_CROSSOVER || deg(v) <= NTL_zz_pX_GCD_CROSSOVER) { 
      PlainResultant(rres, u, v);
      return;
   }

   zz_pX u1, v1;

   u1 = u;
   v1 = v;

   zz_p res, t;
   set(res);

   if (deg(u1) == deg(v1)) {
      rem(u1, u1, v1);
      swap(u1, v1);

      if (IsZero(v1)) {
         clear(rres);
         return;
      }

      power(t, LeadCoeff(u1), deg(u1) - deg(v1));
      mul(res, res, t);
      if (deg(u1) & 1)
         negate(res, res);
   }
   else if (deg(u1) < deg(v1)) {
      swap(u1, v1);
      if (deg(u1) & deg(v1) & 1)
         negate(res, res);
   }

   // deg(u1) > deg(v1) && v1 != 0

   vec_zz_p cvec;
   vec_long  dvec;

   cvec.SetMaxLength(deg(v1)+2);
   dvec.SetMaxLength(deg(v1)+2);

   append(cvec, LeadCoeff(u1));
   append(dvec, deg(u1));


   while (deg(u1) > NTL_zz_pX_GCD_CROSSOVER && !IsZero(v1)) { 
      ResHalfGCD(u1, v1, cvec, dvec);

      if (!IsZero(v1)) {
         append(cvec, LeadCoeff(v1));
         append(dvec, deg(v1));
         rem(u1, u1, v1);
         swap(u1, v1);
      }
   }

   if (IsZero(v1) && deg(u1) > 0) {
      clear(rres);
      return;
   }

   long i, l;
   l = dvec.length();

   if (deg(u1) == 0) {
      // we went all the way...

      for (i = 0; i <= l-3; i++) {
         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
         mul(res, res, t);
         if (dvec[i] & dvec[i+1] & 1)
            negate(res, res);
      }

      power(t, cvec[l-1], dvec[l-2]);
      mul(res, res, t);
   }
   else {
      for (i = 0; i <= l-3; i++) {
         power(t, cvec[i+1], dvec[i]-dvec[i+2]);
         mul(res, res, t);
         if (dvec[i] & dvec[i+1] & 1)
            negate(res, res);
      }

      power(t, cvec[l-1], dvec[l-2]-deg(v1));
      mul(res, res, t);
      if (dvec[l-2] & dvec[l-1] & 1)
         negate(res, res);

      PlainResultant(t, u1, v1);
      mul(res, res, t);
   }

   rres = res;
}

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

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

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

   x = t;
}


NTL_END_IMPL