#include <NTL/ZZ_pX.h>#include <NTL/new.h>NTL_START_IMPLlong 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 (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 (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(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.   // We optimize this case, as it does sometimes arise naturally   // in some situations.   long n = (1L << k);   long xx;   ZZ_p a0, a1, b0, b1, c0, d0, u0, u1, v0, v1, nu0, nu1, nv0;   static ZZ 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, rep(a0), rep(u0));      mul(t2, rep(b0), rep(v0));      add(t1, t1, t2);       conv(nu0, t1);      mul(t1, rep(a1), rep(u0));      mul(t2, rep(a0), rep(u1));      add(t1, t1, t2);      mul(t2, rep(b1), rep(v0));      add(t1, t1, t2);      mul(t2, rep(b0), rep(v1));      add(t1, t1, t2);      conv(nu1, t1);      mul(t1, rep(c0), rep(u0));      mul(t2, rep(d0), rep(v0));      add (t1, t1, t2);      conv(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, rep(a0), rep(u0));      mul(t2, rep(b0), rep(v0));      add(t1, t1, t2);       conv(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;   ZZVec tmp(deg(U)+1, ZZ_pInfo->ExtendedModulusSize);   ZZ_pX Q, t(INIT_SIZE, d_red);   while (deg(V) > goal) {      PlainDivRem(Q, U, U, V, tmp);      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);