/**************************************************************************\MODULE: ZZXSUMMARY:The class ZZX implements polynomials in ZZ[X], i.e., univariatepolynomials with integer coefficients.Polynomial multiplication is implemented using one of 4 differentalgorithms:1) classical 2) Karatsuba3) Schoenhage & Strassen --- performs an FFT by working     modulo a "Fermat number" of appropriate size...     good for polynomials with huge coefficients     and moderate degree4) CRT/FFT --- performs an FFT by working modulo several     small primes...good for polynomials with moderate coefficients     and huge degree.The choice of algorithm is somewhat heuristic, and may not always beperfect.Many thanks to Juergen Gerhard <jngerhar@plato.uni-paderborn.de> forpointing out the deficiency in the NTL-1.0 ZZX arithmetic, and forcontributing the Schoenhage/Strassen code.Extensive use is made of modular algorithms to enhance performance(e.g., the GCD algorithm and amny others).\**************************************************************************/#include <NTL/vec_ZZ.h>#include "zz_pX.h"#include <NTL/ZZ_pX.h>class ZZX {public:   ZZX(); // initial value 0   ZZX(const ZZX& a); // copy   ZZX& operator=(const ZZX& a); // assignment   ZZX& operator=(const ZZ& a);   ZZX& operator=(long a);   ~ZZX();   ZZX(long i, const ZZ& c); // initial value X^i*c   ZZX(long i, long c);    // ...};/**************************************************************************\                                  Comparison\**************************************************************************/long operator==(const ZZX& a, const ZZX& b);long operator!=(const ZZX& a, const ZZX& b);long IsZero(const ZZX& a);  // test for 0long IsOne(const ZZX& a);  // test for 1// PROMOTIONS: operators ==, != promote {long, ZZ} to ZZX on (a, b)./**************************************************************************\                                   Addition\**************************************************************************/// operator notation:ZZX operator+(const ZZX& a, const ZZX& b);ZZX operator-(const ZZX& a, const ZZX& b);ZZX operator-(const ZZX& a); // unary -ZZX& operator+=(ZZX& x, const ZZX& a);ZZX& operator-=(ZZX& x, const ZZX& a);ZZX& operator++(ZZX& x);  // prefixvoid operator++(ZZX& x, int);  // postfixZZX& operator--(ZZX& x);  // prefixvoid operator--(ZZX& x, int);  // postfix// procedural versions:void add(ZZX& x, const ZZX& a, const ZZX& b); // x = a + bvoid sub(ZZX& x, const ZZX& a, const ZZX& b); // x = a - bvoid negate(ZZX& x, const ZZX& a); // x = -a// PROMOTIONS: binary +, - and procedures add, sub promote {long, ZZ} // to ZZX on (a, b)./**************************************************************************\                               Multiplication\**************************************************************************/// operator notation:ZZX operator*(const ZZX& a, const ZZX& b);ZZX& operator*=(ZZX& x, const ZZX& a);// procedural versions:void mul(ZZX& x, const ZZX& a, const ZZX& b); // x = a * bvoid sqr(ZZX& x, const ZZX& a); // x = a^2ZZX sqr(const ZZX& a);// PROMOTIONS: operator * and procedure mul promote {long, ZZ} to ZZX // on (a, b)./**************************************************************************\                               Shift OperationsLeftShift by n means multiplication by X^nRightShift by n means division by X^nA negative shift amount reverses the direction of the shift.\**************************************************************************/// operator notation:ZZX operator<<(const ZZX& a, long n);ZZX operator>>(const ZZX& a, long n);ZZX& operator<<=(ZZX& x, long n);ZZX& operator>>=(ZZX& x, long n);// procedural versions:void LeftShift(ZZX& x, const ZZX& a, long n); ZZX LeftShift(const ZZX& a, long n);void RightShift(ZZX& x, const ZZX& a, long n); ZZX RightShift(const ZZX& a, long n); /**************************************************************************\                                  Division\**************************************************************************/// Given polynomials a, b in ZZ[X], there exist polynomials// q, r in QQ[X] such that a = b*q + r, deg(r) < deg(b).// These routines return q and/or r if q and/or r lie(s) in ZZ[X],// and otherwise raise an error.  // Note that if the leading coefficient of b is 1 or -1, // then q and r always lie in ZZ[X], and no error can occur.// For example, you can write f/2 for a ZZX f.  If all coefficients// of f are even, the result is f with a factor of two removed;// otherwise, an error is raised.  More generally, f/g will be// evaluate q in ZZ[X] such that f = q*g if such a q exists,// and will otherwise raise an error.// See also below the routines for pseudo-division and division// predicates for routines that are perhaps more useful in// some situations.// operator notation: ZZX operator/(const ZZX& a, const ZZX& b);ZZX operator/(const ZZX& a, const ZZ& b);ZZX operator/(const ZZX& a, long b);ZZX operator%(const ZZX& a, const ZZX& b);ZZX& operator/=(ZZX& x, const ZZX& b);ZZX& operator/=(ZZX& x, const ZZ& b);ZZX& operator/=(ZZX& x, long b);ZZX& operator%=(ZZX& x, const ZZX& b);// procedural versions:void DivRem(ZZX& q, ZZX& r, const ZZX& a, const ZZX& b);// computes q, r such that a = b q + r and deg(r) < deg(b).// requires LeadCoeff(b) is a unit (+1, -1); otherwise,// an error is raised.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);// same as DivRem, but only computes qvoid rem(ZZX& r, const ZZX& a, const ZZX& b);// same as DivRem, but only computes r// divide predicates:long divide(ZZX& q, const ZZX& a, const ZZX& b);long divide(ZZX& q, const ZZX& a, const ZZ& b);long divide(ZZX& q, const ZZX& a, long b);// if b | a, sets q = a/b and returns 1; otherwise returns 0long divide(const ZZX& a, const ZZX& b);long divide(const ZZX& a, const ZZ& b);long divide(const ZZX& a, long b);// if b | a, returns 1; otherwise returns 0// These algorithms employ a modular approach, performing the division// modulo small primes (reconstructing q via the CRT).  It is// usually much faster than the general division routines above// (especially when b does not divide a).void content(ZZ& d, const ZZX& f);ZZ content(const ZZX& f);// d = content of f, sign(d) == sign(LeadCoeff(f)); content(0) == 0void PrimitivePart(ZZX& pp, const ZZX& f);ZZX PrimitivePart(const ZZX& f); // pp = primitive part of f, LeadCoeff(pp) >= 0; PrimitivePart(0) == 0// pseudo-division: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.  Only the classical// algorithm is used.void PseudoDiv(ZZX& q, const ZZX& a, const ZZX& b);ZZX PseudoDiv(const ZZX& a, const ZZX& b);// same as PseudoDivRem, but only computes qvoid PseudoRem(ZZX& r, const ZZX& a, const ZZX& b);ZZX PseudoRem(const ZZX& a, const ZZX& b);// same as PseudoDivRem, but only computes r/**************************************************************************\                                  GCD's