#ifndef POLY_H
#define POLY_H 1

#include <assert.h>
#include <iostream>
#include <vector>

/*
  this keeps a rational polynomial as the quotient of an integer polynomial
  and an integer, for efficiency reasons.
  the fraction is not automatically reduced, to optimize addition.
  each poly has a maxdeg (relating to space allocation) and a deg; the deg 
  is -1 for the 0 poly and 0 for other constant polys.
  when a poly is declared its maxdeg is specified; this can be changed 
  using set_maxdeg()

  all coefficients above the degree are kept at 0
*/

class poly {
private:
    vector <INT> numer;
    INT denom;
    int d;
    int maxd;
    void setdeg(int startd) {
	for(d=startd;d>=0 && numer[d] == 0; d--);
    }
public:
    explicit poly(int n) :numer(n+1), denom(1) {d = -1; maxd = n;}
    poly() :numer(1), denom(1) {d = -1; maxd = 0;}
    ~poly() {}
    int maxdeg() const {return maxd;}
    int deg() const {return d;}

// this is expensive, so use sparingly    
    void set_maxdeg(int n);

    bool zero() const {return deg() == -1;}

// operator[] caused too many problems with maintaining the degree 
//  (when used as an lvalue) and was misleading (since it only referred
//  to the numerator

    const INT& get_numer(int k) const {
	assert(0<=k && k<= maxdeg());
	return numer[k];
    }
    const INT& get_denom() const {return denom;}
    void set_denom(const INT& C) {
	assert(C!=0);
	denom = C;
    }
    void set_numer(int k, const INT& C) {
	if(k>deg()) {
	    if(C==0)
		return;
	    assert(k<=maxdeg());
	    numer[k] = C;
	    d = k;
	    return;
	}
	numer[k] = C;
	if(k==deg() && C == 0)
	    setdeg(k);
    }
    void set_numer(int k, int c) {
	set_numer(k,INT(c));
    }


//  unlike the ops in algebra.h, operations do not have unexpected 
//  side effects on the operands

    poly& operator=(int z) {
	for(int k=deg(); k>0; k--)
	    numer[k] = 0;
	numer[0] = z;
	denom = 1;
	d = z == 0 ? -1 : 0;
	return *this;
    }
    poly& operator=(const INT z) {
	for(int k=deg(); k>0; k--)
	    numer[k] = 0;
	numer[0] = z;
	denom = 1;
	d = z == 0 ? -1 : 0;
	return *this;
    }
    poly& operator=(const poly& Q) {
	assert(maxdeg()>=Q.deg());
	denom = Q.denom;
	for(int k=Q.d;k>=0;k--)
	    numer[k] = Q.numer[k];
	for(int k=Q.d+1;k<=deg();k++)
	    numer[k] = 0;
	d = Q.d;
	return *this;
    }
    poly& operator*=(int z) {
	if(z == 0)
	    *this = 0;
	else 
	    for(int k = deg(); k>=0; k--)
		numer[k] *= z;
	return *this;
    }
    poly& operator*=(const INT& z) {
	if(z == 0)
	    *this = 0;
	else
	    for(int k = deg(); k>=0; k--)
		numer[k] *= z;
	return *this;
    }	
    poly& operator*=(const poly& Q) {
	if (zero() || Q.zero())
	    return *this = 0;
	int newd = d + Q.d;
	assert(newd <= maxdeg());
	denom *= Q.denom;
	INT acc;
	for(int k=newd;k>=0;k--) {
	    acc = 0;
	    int starti = k - Q.d;
	    if (starti < 0)
		starti = 0;
	    int endi = k < d ? k : d;
	    for(int i=endi; i>= starti; i--)
		acc += numer[i]*Q.numer[k-i];
	    numer[k] = acc;
	}
	d = newd;
	return *this;
    }
/*      poly& pow(poly& V, uint p) const { */
/*  	V = *this; */
/*  	for(int k = 2; k<=p; k++) */
/*  	    V *= *this; */
/*  	return V; */
/*      } */
    poly& operator^= (int p) {
	assert(p >= 0);
	if(zero() || p == 1)
	    return *this;
	if(p == 0)
	    return *this = 1;
	poly V(d);
	V = *this;
	for(int k = 2; k<=p; k++)
	    *this *= V;
	return *this;
    }
    poly& operator/=(int z) {
	assert(z!=0);
	denom *= z;
	return *this;
    }	
    poly& operator/=(const INT& z) {
	assert(z!=0);
	denom *= z;
	return *this;
    }	
    poly& operator+=(const poly& Q) {
	assert(maxdeg() >= Q.deg());
	if(denom == Q.denom)
	    for(int k=Q.deg();k>=0;k--)
		numer[k] += Q.numer[k];
	else {
	    INT GCD;
	    gcd(GCD,denom,Q.denom);
	    INT x,y;
	    x = Q.denom/GCD;
	    y = denom/GCD;
	    *this *= x;
	    denom *= x;
	    for(int k=Q.deg();k>=0;k--)
		numer[k] += y*Q.numer[k];
	}
	    
	if (Q.deg() > deg())
	    d = Q.deg();
	else if (Q.deg() == deg() && numer[d] == 0)
	    setdeg(d);
	return *this;
    }
    RAT& eval(RAT& y,int x) const {
	INT z = 0;
	for(int k=deg(); k>= 0; k--)
	    z = z*x + numer[k];
	fraction(y,z,denom);
	return y;
    }
    
// cancels common factors between numer and denom
    poly& reduce();

// arranges that denom == D
// returns false if not possible (i.e., if D is not a multiple of the 
//   reduced denom
    bool fixdenom(const INT& D);
    bool fixdenom(int d) {
	return fixdenom(INT(d));
    }

    friend ostream& operator<<(ostream& f, const poly& z);
    friend istream& operator>>(istream& f, poly& z);
};

#endif