inline ZZX operator-(const ZZX& a, const ZZ& b)
   { ZZX x; sub(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator-(const ZZX& a, long b)
   { ZZX x; sub(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator-(const ZZ& a, const ZZX& b)
   { ZZX x; sub(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator-(long a, const ZZX& b)
   { ZZX x; sub(x, a, b); NTL_OPT_RETURN(ZZX, x); }


inline ZZX& operator+=(ZZX& x, const ZZX& b)
   { add(x, x, b); return x; }

inline ZZX& operator+=(ZZX& x, const ZZ& b)
   { add(x, x, b); return x; }

inline ZZX& operator+=(ZZX& x, long b)
   { add(x, x, b); return x; }

inline ZZX& operator-=(ZZX& x, const ZZX& b)
   { sub(x, x, b); return x; }

inline ZZX& operator-=(ZZX& x, const ZZ& b)
   { sub(x, x, b); return x; }

inline ZZX& operator-=(ZZX& x, long b)
   { sub(x, x, b); return x; }


inline ZZX operator-(const ZZX& a) 
   { ZZX x; negate(x, a); NTL_OPT_RETURN(ZZX, x); }

inline ZZX& operator++(ZZX& x) { add(x, x, 1); return x; }
inline void operator++(ZZX& x, int) { add(x, x, 1); }
inline ZZX& operator--(ZZX& x) { sub(x, x, 1); return x; }
inline void operator--(ZZX& x, int) { sub(x, x, 1); }


/*****************************************************************

                        Multiplication

******************************************************************/


void mul(ZZX& x, const ZZX& a, const ZZX& b);
// x = a * b


void sqr(ZZX& x, const ZZX& a);
inline ZZX sqr(const ZZX& a)
   { ZZX x; sqr(x, a); NTL_OPT_RETURN(ZZX, x); }
// x = a^2

void PlainMul(ZZX& x, const ZZX& a, const ZZX& b);
void PlainSqr(ZZX& x, const ZZX& a);

void KarMul(ZZX& x, const ZZX& a, const ZZX& b);
void KarSqr(ZZX& x, const ZZX& a);

void HomMul(ZZX& x, const ZZX& a, const ZZX& b);
void HomSqr(ZZX& x, const ZZX& a);

void SSMul(ZZX& x, const ZZX& a, const ZZX& b);
void SSSqr(ZZX& x, const ZZX& a);

double SSRatio(long na, long maxa, long nb, long maxb);


void mul(ZZX & x, const ZZX& a, const ZZ& b);
void mul(ZZX& x, const ZZX& a, long b);

inline void mul(ZZX& x, const ZZ& a, const ZZX& b) { mul(x, b, a); } 
inline void mul(ZZX& x, long a, const ZZX& b) { mul(x, b, a); } 


inline ZZX operator*(const ZZX& a, const ZZX& b)
   { ZZX x; mul(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator*(const ZZX& a, const ZZ& b)
   { ZZX x; mul(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator*(const ZZX& a, long b)
   { ZZX x; mul(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator*(const ZZ& a, const ZZX& b)
   { ZZX x; mul(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator*(long a, const ZZX& b)
   { ZZX x; mul(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX& operator*=(ZZX& x, const ZZX& b)
   { mul(x, x, b); return x; }

inline ZZX& operator*=(ZZX& x, const ZZ& b)
   { mul(x, x, b); return x; }

inline ZZX& operator*=(ZZX& x, long b)
   { mul(x, x, b); return x; }






/*************************************************************

                      Division

**************************************************************/



// "plain" versions
void PlainPseudoDivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b);
void PlainPseudoDiv(ZZX& q, const ZZX& a, const ZZX& b);
void PlainPseudoRem(ZZX& r, const ZZX& a, const ZZX& b);

// "homomorphic imaging" versions
void HomPseudoDivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b);
void HomPseudoDiv(ZZX& q, const ZZX& a, const ZZX& b);
void HomPseudoRem(ZZX& r, const ZZX& a, const ZZX& b);

inline void PseudoDivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b)
// performs pseudo-division:  computes q and r
// with deg(r) < deg(b), and LeadCoeff(b)^(deg(a)-deg(b)+1) a = b q + r.
// current implementation always defaults to "plain"

   { PlainPseudoDivRem(q, r, a, b); }

inline void PseudoDiv(ZZX& q, const ZZX& a, const ZZX& b)

   { PlainPseudoDiv(q, a, b); }

inline void PseudoRem(ZZX& r, const ZZX& a, const ZZX& b)

   { PlainPseudoRem(r, a, b); }

inline ZZX PseudoDiv(const ZZX& a, const ZZX& b)
   { ZZX x; PseudoDiv(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX PseudoRem(const ZZX& a, const ZZX& b)
   { ZZX x; PseudoRem(x, a, b); NTL_OPT_RETURN(ZZX, x); }


#ifndef NTL_TRANSITION

void DivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b);

void div(ZZX& q, const ZZX& a, const ZZX& b);
void div(ZZX& q, const ZZX& a, const ZZ& b);
void div(ZZX& q, const ZZX& a, long b);

void rem(ZZX& r, const ZZX& a, const ZZX& b);

inline ZZX operator/(const ZZX& a, const ZZX& b)
   { ZZX x; div(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator/(const ZZX& a, const ZZ& b)
   { ZZX x; div(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX operator/(const ZZX& a, long b)
   { ZZX x; div(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX& operator/=(ZZX& x, const ZZ& b)
   { div(x, x, b); return x; }

inline ZZX& operator/=(ZZX& x, long b)
   { div(x, x, b); return x; }

inline ZZX& operator/=(ZZX& x, const ZZX& b)
   { div(x, x, b); return x; }


inline ZZX operator%(const ZZX& a, const ZZX& b)
   { ZZX x; rem(x, a, b); NTL_OPT_RETURN(ZZX, x); }

inline ZZX& operator%=(ZZX& x, const ZZX& b)
   { rem(x, x, b); return x; }

#endif


// Modular arithemtic---f must be monic, and other args
// must have degree less than that of f

void MulMod(ZZX& x, const ZZX& a, const ZZX& b, const ZZX& f);

inline ZZX MulMod(const ZZX& a, const ZZX& b, const ZZX& f)
   { ZZX x; MulMod(x, a, b, f); NTL_OPT_RETURN(ZZX, x); }

void SqrMod(ZZX& x, const ZZX& a, const ZZX& f);

inline ZZX SqrMod(const ZZX& a, const ZZX& f)
   { ZZX x; SqrMod(x, a, f); NTL_OPT_RETURN(ZZX, x); }

void MulByXMod(ZZX& x, const ZZX& a, const ZZX& f);

inline ZZX MulByXMod(const ZZX& a, const ZZX& f)
   { ZZX x; MulByXMod(x, a, f); NTL_OPT_RETURN(ZZX, x); }


// these always use "plain" division
long PlainDivide(ZZX& q, const ZZX& a, const ZZX& b);
long PlainDivide(const ZZX& a, const ZZX& b);

// these always use "homomorphic imaging"
long HomDivide(ZZX& q, const ZZX& a, const ZZX& b);
long HomDivide(const ZZX& a, const ZZX& b);

long divide(ZZX& q, const ZZX& a, const ZZX& b);
// if b | a, sets q = a/b and returns 1; otherwise returns 0

long divide(const ZZX& a, const ZZX& b);


long divide(ZZX& q, const ZZX& a, const ZZ& b);
// if b | a, sets q = a/b and returns 1; otherwise returns 0

long divide(const ZZX& a, const ZZ& b);
// if b | a, returns 1; otherwise returns 0

//single-precision versions
long divide(ZZX& q, const ZZX& a, long b);
long divide(const ZZX& a, long b);



void content(ZZ& d, const ZZX& f);
// d = content of f, sign(d) == sign(LeadCoeff(f))

inline ZZ content(const ZZX& f)
   { ZZ x; content(x, f); NTL_OPT_RETURN(ZZ, x); }



void PrimitivePart(ZZX& pp, const ZZX& f);
// pp = primitive part of f, LeadCoeff(pp) >= 0

inline ZZX PrimitivePart(const ZZX& f)
   { ZZX x; PrimitivePart(x, f); NTL_OPT_RETURN(ZZX, x); }
   

void GCD(ZZX& d, const ZZX& a, const ZZX& b);
// d = gcd(a, b), LeadCoeff(d) >= 0

inline ZZX GCD(const ZZX& a, const ZZX& b)
   { ZZX x; GCD(x, a, b); NTL_OPT_RETURN(ZZX, x); }

long MaxBits(const ZZX& f);
// returns max NumBits of coefficients of f

long CharPolyBound(const ZZX& a, const ZZX& f);



/***************************************************************

                      traces, norms, resultants

****************************************************************/

void TraceVec(vec_ZZ& S, const ZZX& f);
// S[i] = Trace(X^i mod f), for i = 0..deg(f)-1.
// f must be a monic polynomial.

inline vec_ZZ TraceVec(const ZZX& f)
   { vec_ZZ x; TraceVec(x, f); NTL_OPT_RETURN(vec_ZZ, x); }

void TraceMod(ZZ& res, const ZZX& a, const ZZX& f);
inline ZZ TraceMod(const ZZX& a, const ZZX& f)
   { ZZ x; TraceMod(x, a, f); NTL_OPT_RETURN(ZZ, x); }
// res = trace of (a mod f)
// f must be monic


void resultant(ZZ& res, const ZZX& a, const ZZX& b, long deterministic=0);
inline ZZ resultant(const ZZX& a, const ZZX& b, long deterministic=0)
   { ZZ x; resultant(x, a, b, deterministic); NTL_OPT_RETURN(ZZ, x); }

// res = resultant of a and b
// if !deterministic, then it may use a randomized strategy
//    that errs with probability no more than 2^{-80}.

void NormMod(ZZ& res, const ZZX& a, const ZZX& f, long deterministic=0);
inline ZZ NormMod(const ZZX& a, const ZZX& f, long deterministic=0)
   { ZZ x; NormMod(x, a, f, deterministic); NTL_OPT_RETURN(ZZ, x); }
// res = norm of (a mod f)
// f must be monic
// if !deterministic, then it may use a randomized strategy
//    that errs with probability no more than 2^{-80}.


void discriminant(ZZ& d, const ZZX& a, long deterministic=0);
inline ZZ discriminant(const ZZX& a, long deterministic=0)
   { ZZ x; discriminant(x, a, deterministic); NTL_OPT_RETURN(ZZ, x); }
// d = discriminant of a
//   = (-1)^{m(m-1)/2} resultant(a, a')/lc(a),
//     where m = deg(a)
// if !deterministic, then it may use a randomized strategy
//    that errs with probability no more than 2^{-80}.


void CharPolyMod(ZZX& g, const ZZX& a, const ZZX& f, long deterministic=0);
inline ZZX CharPolyMod(const ZZX& a, const ZZX& f, long deterministic=0)
   { ZZX x; CharPolyMod(x, a, f, deterministic); NTL_OPT_RETURN(ZZX, x); }
// g = char poly of (a mod f)
// f must be monic
// if !deterministic, then it may use a randomized strategy
//    that errs with probability no more than 2^{-80}.


void MinPolyMod(ZZX& g, const ZZX& a, const ZZX& f);
inline ZZX MinPolyMod(const ZZX& a, const ZZX& f)
   { ZZX x; MinPolyMod(x, a, f); NTL_OPT_RETURN(ZZX, x); }
// g = min poly of (a mod f)
// f must be monic
// may use a probabilistic strategy that errs with
//   probability no more than 2^{-80}


void XGCD(ZZ& r, ZZX& s, ZZX& t, const ZZX& a, const ZZX& b, 
          long deterministic=0);
// r = resultant of a and b;
// if r != 0, then computes s and t such that:
//   a*s + b*t = r;
// otherwise s and t not affected.
// if !deterministic, then resultant computation may use a randomized strategy
//    that errs with probability no more than 2^{-80}.



/******************************************************

      Incremental Chinese Remaindering

*******************************************************/

long CRT(ZZX& a, ZZ& prod, const zz_pX& A);
long CRT(ZZX& a, ZZ& prod, const ZZ_pX& A);
// If p is the current modulus with (p, prod) = 1;
// Computes b such that b = a mod prod and b = A mod p,
//    with coefficients in the interval (-p*prod/2, p*prod/2];
// Sets a = b, prod = p*prod, and returns 1 if a's value changed.




// vectors

NTL_vector_decl(ZZX,vec_ZZX)

NTL_eq_vector_decl(ZZX,vec_ZZX)

NTL_io_vector_decl(ZZX,vec_ZZX)

NTL_CLOSE_NNS

#endif