/**************************************************************************\

MODULE: ZZX

SUMMARY:

The class ZZX implements polynomials in ZZ[X], i.e., univariate
polynomials with integer coefficients.

Polynomial multiplication is implemented using one of 4 different
algorithms:

1) classical 

2) Karatsuba

3) Schoenhage & Strassen --- performs an FFT by working
     modulo a "Fermat number" of appropriate size...
     good for polynomials with huge coefficients
     and moderate degree

4) 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 be
perfect.

Many thanks to Juergen Gerhard <jngerhar@plato.uni-paderborn.de> for
pointing out the deficiency in the NTL-1.0 ZZX arithmetic, and for
contributing 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 0
long 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);  // prefix
void operator++(ZZX& x, int);  // postfix

ZZX& operator--(ZZX& x);  // prefix
void operator--(ZZX& x, int);  // postfix


// procedural versions:

void add(ZZX& x, const ZZX& a, const ZZX& b); // x = a + b
void sub(ZZX& x, const ZZX& a, const ZZX& b); // x = a - b
void 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 * b

void sqr(ZZX& x, const ZZX& a); // x = a^2
ZZX sqr(const ZZX& a);

// PROMOTIONS: operator * and procedure mul promote {long, ZZ} to ZZX 
// on (a, b).


/**************************************************************************\

                               Shift Operations

LeftShift by n means multiplication by X^n
RightShift by n means division by X^n

A 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 q

void 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 0


long 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) == 0

void 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 q

void 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