/*****************************************************************                        Multiplication******************************************************************/void mul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b);// x = a * bvoid sqr(ZZ_pX& x, const ZZ_pX& a);inline ZZ_pX sqr(const ZZ_pX& a)   { ZZ_pX x; sqr(x, a); NTL_OPT_RETURN(ZZ_pX, x); }// x = a^2void mul(ZZ_pX & x, const ZZ_pX& a, const ZZ_p& b); void mul(ZZ_pX& x, const ZZ_pX& a, long b);inline void mul(ZZ_pX& x, const ZZ_p& a, const ZZ_pX& b)   { mul(x, b, a); }inline void mul(ZZ_pX& x, long a, const ZZ_pX& b)   { mul(x, b, a); }inline ZZ_pX operator*(const ZZ_pX& a, const ZZ_pX& b)   { ZZ_pX x; mul(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX operator*(const ZZ_pX& a, const ZZ_p& b)   { ZZ_pX x; mul(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX operator*(const ZZ_pX& a, long b)   { ZZ_pX x; mul(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX operator*(const ZZ_p& a, const ZZ_pX& b)   { ZZ_pX x; mul(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX operator*(long a, const ZZ_pX& b)   { ZZ_pX x; mul(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX& operator*=(ZZ_pX& x, const ZZ_pX& b)   { mul(x, x, b); return x; }inline ZZ_pX& operator*=(ZZ_pX& x, const ZZ_p& b)   { mul(x, x, b); return x; }inline ZZ_pX& operator*=(ZZ_pX& x, long b)   { mul(x, x, b); return x; }void PlainMul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b);// always uses the "classical" algorithmvoid PlainSqr(ZZ_pX& x, const ZZ_pX& a);// always uses the "classical" algorithmvoid FFTMul(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b);// always uses the FFTvoid FFTSqr(ZZ_pX& x, const ZZ_pX& a);// always uses the FFTvoid MulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n);// x = a * b % X^ninline ZZ_pX MulTrunc(const ZZ_pX& a, const ZZ_pX& b, long n)   { ZZ_pX x; MulTrunc(x, a, b, n); NTL_OPT_RETURN(ZZ_pX, x); }void PlainMulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n);void FFTMulTrunc(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, long n);void SqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n);// x = a^2 % X^ninline ZZ_pX SqrTrunc(const ZZ_pX& a, long n)   { ZZ_pX x; SqrTrunc(x, a, n); NTL_OPT_RETURN(ZZ_pX, x); }void PlainSqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n);void FFTSqrTrunc(ZZ_pX& x, const ZZ_pX& a, long n);void power(ZZ_pX& x, const ZZ_pX& a, long e);inline ZZ_pX power(const ZZ_pX& a, long e)   { ZZ_pX x; power(x, a, e); NTL_OPT_RETURN(ZZ_pX, x); }// The following data structures and routines allow one// to hand-craft various algorithms, using the FFT convolution// algorithms directly.// Look in the file ZZ_pX.c for examples.// FFT representation of polynomialsclass FFTRep {public:   long k;                // a 2^k point representation   long MaxK;             // maximum space allocated   long **tbl;   long NumPrimes;    void SetSize(long NewK);   FFTRep(const FFTRep& R);   FFTRep& operator=(const FFTRep& R);   FFTRep() { k = MaxK = -1; tbl = 0; NumPrimes = 0; }   FFTRep(INIT_SIZE_TYPE, long InitK)    { k = MaxK = -1; tbl = 0; NumPrimes = 0; SetSize(InitK); }   ~FFTRep();};void ToFFTRep(FFTRep& y, const ZZ_pX& x, long k, long lo, long hi);// computes an n = 2^k point convolution of x[lo..hi].inline void ToFFTRep(FFTRep& y, const ZZ_pX& x, long k)   { ToFFTRep(y, x, k, 0, deg(x)); }void RevToFFTRep(FFTRep& y, const vec_ZZ_p& x,                 long k, long lo, long hi, long offset);// computes an n = 2^k point convolution of X^offset*x[lo..hi] mod X^n-1// using "inverted" evaluation points.void FromFFTRep(ZZ_pX& x, FFTRep& y, long lo, long hi);// converts from FFT-representation to coefficient representation// only the coefficients lo..hi are computed// NOTE: this version destroys the data in y// non-destructive versions of the abovevoid NDFromFFTRep(ZZ_pX& x, const FFTRep& y, long lo, long hi, FFTRep& temp);void NDFromFFTRep(ZZ_pX& x, const FFTRep& y, long lo, long hi);void RevFromFFTRep(vec_ZZ_p& x, FFTRep& y, long lo, long hi);   // converts from FFT-representation to coefficient representation   // using "inverted" evaluation points.   // only the coefficients lo..hi are computedvoid FromFFTRep(ZZ_p* x, FFTRep& y, long lo, long hi);// convert out coefficients lo..hi of y, store result in x.// no normalization is done.// direct manipulation of FFT repsvoid mul(FFTRep& z, const FFTRep& x, const FFTRep& y);void sub(FFTRep& z, const FFTRep& x, const FFTRep& y);void add(FFTRep& z, const FFTRep& x, const FFTRep& y);void reduce(FFTRep& x, const FFTRep& a, long k);// reduces a 2^l point FFT-rep to a 2^k point FFT-repvoid AddExpand(FFTRep& x, const FFTRep& a);//  x = x + (an "expanded" version of a)// This data structure holds unconvoluted modular representations// of polynomialsclass ZZ_pXModRep {private:   ZZ_pXModRep(const ZZ_pXModRep&); // disabled   void operator=(const ZZ_pXModRep&); // disabledpublic:   long n;   long MaxN;   long **tbl;   long NumPrimes;    void SetSize(long NewN);   ZZ_pXModRep() { n = MaxN = 0; tbl = 0; NumPrimes = 0; }   ZZ_pXModRep(INIT_SIZE_TYPE, long k)    { n = MaxN = 0; tbl = 0; NumPrimes = 0; SetSize(k); }   ~ZZ_pXModRep();};void ToZZ_pXModRep(ZZ_pXModRep& x, const ZZ_pX& a, long lo, long hi);void ToFFTRep(FFTRep& x, const ZZ_pXModRep& a, long k, long lo, long hi);// converts coefficients lo..hi to a 2^k-point FFTRep.// must have hi-lo+1 < 2^k/*************************************************************                      Division**************************************************************/void DivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b);// q = a/b, r = a%bvoid div(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b);// q = a/bvoid div(ZZ_pX& q, const ZZ_pX& a, const ZZ_p& b);void div(ZZ_pX& q, const ZZ_pX& a, long b);void rem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b);// r = a%blong divide(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b);// if b | a, sets q = a/b and returns 1; otherwise returns 0long divide(const ZZ_pX& a, const ZZ_pX& b);// if b | a, sets q = a/b and returns 1; otherwise returns 0void InvTrunc(ZZ_pX& x, const ZZ_pX& a, long m);// computes x = a^{-1} % X^m // constant term must be non-zeroinline ZZ_pX InvTrunc(const ZZ_pX& a, long m)   { ZZ_pX x; InvTrunc(x, a, m); NTL_OPT_RETURN(ZZ_pX, x); }// These always use "classical" arithmeticvoid PlainDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b);void PlainDiv(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b);void PlainRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b);void PlainRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b, ZZVec& tmp);void PlainDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b,                 ZZVec& tmp);// These always use FFT arithmeticvoid FFTDivRem(ZZ_pX& q, ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b);void FFTDiv(ZZ_pX& q, const ZZ_pX& a, const ZZ_pX& b);void FFTRem(ZZ_pX& r, const ZZ_pX& a, const ZZ_pX& b);void PlainInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m);// always uses "classical" algorithm// ALIAS RESTRICTION: input may not alias outputvoid NewtonInvTrunc(ZZ_pX& x, const ZZ_pX& a, long m);// uses a Newton Iteration with the FFT.// ALIAS RESTRICTION: input may not alias outputinline ZZ_pX operator/(const ZZ_pX& a, const ZZ_pX& b)   { ZZ_pX x; div(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX operator/(const ZZ_pX& a, const ZZ_p& b)   { ZZ_pX x; div(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX operator/(const ZZ_pX& a, long b)   { ZZ_pX x; div(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX& operator/=(ZZ_pX& x, const ZZ_p& b)   { div(x, x, b); return x; }inline ZZ_pX& operator/=(ZZ_pX& x, long b)   { div(x, x, b); return x; }inline ZZ_pX& operator/=(ZZ_pX& x, const ZZ_pX& b)   { div(x, x, b); return x; }inline ZZ_pX operator%(const ZZ_pX& a, const ZZ_pX& b)   { ZZ_pX x; rem(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }inline ZZ_pX& operator%=(ZZ_pX& x, const ZZ_pX& b)   { rem(x, x, b); return x; }/***********************************************************                         GCD's************************************************************/void GCD(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b);// x = GCD(a, b),  x is always monic (or zero if a==b==0).inline ZZ_pX GCD(const ZZ_pX& a, const ZZ_pX& b)   { ZZ_pX x; GCD(x, a, b); NTL_OPT_RETURN(ZZ_pX, x); }void XGCD(ZZ_pX& d, ZZ_pX& s, ZZ_pX& t, const ZZ_pX& a, const ZZ_pX& b);// d = gcd(a,b), a s + b t = d void PlainXGCD(ZZ_pX& d, ZZ_pX& s, ZZ_pX& t, const ZZ_pX& a, const ZZ_pX& b);// same as above, but uses classical algorithmvoid PlainGCD(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b);// always uses "cdlassical" arithmeticclass ZZ_pXMatrix {private:   ZZ_pXMatrix(const ZZ_pXMatrix&);  // disable   ZZ_pX elts[2][2];public:   ZZ_pXMatrix() { }   ~ZZ_pXMatrix() { }   void operator=(const ZZ_pXMatrix&);   ZZ_pX& operator() (long i, long j) { return elts[i][j]; }   const ZZ_pX& operator() (long i, long j) const { return elts[i][j]; }};void HalfGCD(ZZ_pXMatrix& M_out, const ZZ_pX& U, const ZZ_pX& V, long d_red);// deg(U) > deg(V),   1 <= d_red <= deg(U)+1.//// This computes a 2 x 2 polynomial matrix M_out such that//    M_out * (U, V)^T = (U', V')^T,// where U', V' are consecutive polynomials in the Euclidean remainder// sequence of U, V, and V' is the polynomial of highest degree// satisfying deg(V') <= deg(U) - d_red.void XHalfGCD(ZZ_pXMatrix& M_out, ZZ_pX& U, ZZ_pX& V, long d_red);// same as above, except that U is replaced by U', and V by V'/*************************************************************             Modular Arithmetic without pre-conditioning**************************************************************/// arithmetic mod f.// all inputs and outputs are polynomials of degree less than deg(f).// ASSUMPTION: f is assumed monic, and deg(f) > 0.// NOTE: if you want to do many computations with a fixed f,//       use the ZZ_pXModulus data structure and associated routines below.void MulMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& b, const ZZ_pX& f);// x = (a * b) % finline ZZ_pX MulMod(const ZZ_pX& a, const ZZ_pX& b, const ZZ_pX& f)   { ZZ_pX x; MulMod(x, a, b, f); NTL_OPT_RETURN(ZZ_pX, x); }void SqrMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f);// x = a^2 % finline ZZ_pX SqrMod(const ZZ_pX& a,  const ZZ_pX& f)   { ZZ_pX x; SqrMod(x, a, f); NTL_OPT_RETURN(ZZ_pX, x); }void MulByXMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f);// x = (a * X) mod finline ZZ_pX MulByXMod(const ZZ_pX& a,  const ZZ_pX& f)   { ZZ_pX x; MulByXMod(x, a, f); NTL_OPT_RETURN(ZZ_pX, x); }void InvMod(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f);// x = a^{-1} % f, error is a is not invertibleinline ZZ_pX InvMod(const ZZ_pX& a,  const ZZ_pX& f)   { ZZ_pX x; InvMod(x, a, f); NTL_OPT_RETURN(ZZ_pX, x); }long InvModStatus(ZZ_pX& x, const ZZ_pX& a, const ZZ_pX& f);// if (a, f) = 1, returns 0 and sets x = a^{-1} % f// otherwise, returns 1 and sets x = (a, f)/******************************************************************        Modular Arithmetic with Pre-conditioning*******************************************************************/// If you need to do a lot of arithmetic modulo a fixed f,// build ZZ_pXModulus F for f.  This pre-computes information about f// that speeds up the computation a great deal.class ZZ_pXModulus {public:   ZZ_pXModulus() : UseFFT(0), n(-1)  { }   ~ZZ_pXModulus() { }      // the following members may become private in future   ZZ_pX f;   // the modulus   long UseFFT;// flag indicating whether FFT should be used.   long n;     // n = deg(f)   long k;     // least k s/t 2^k >= n   long l;     // least l s/t 2^l >= 2n-3   FFTRep FRep; // 2^k point rep of f                // H = rev((rev(f))^{-1} rem X^{n-1})   FFTRep HRep; // 2^l point rep of H