#ifndef ALGEBRA_H
#define ALGEBRA_H 1

#include <assert.h>
#include <iostream>
#include <sstream>
#include <list>
#include <string>

using namespace std;

class generator;
class expr;

#define ODD(x) (x&1)

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

NOTE: Many of the binary operations destroy their RHS 
(unless of a simple type, like int).
The only RHS variables that are preserved are those declared const.

The main exceptions are =, which always copies, and comparison operators.

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

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

integrands are built from factorlists.
the template expects that class T comes with 
  comparison operators (just < and ==) used for sorting the list of factors
  a unit() function, to indicate which objects can be considered as 1
  a *= operator, which only needs to work for objects that compare ==
  a ^= operator
  an exp() function

objects are considered identical if they test == and have the same 
  exp() values; this is just in the == operator, which is not currently used

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

template <class T> class factorlist: public list<T> {
public:
    factorlist<T>& operator*=(const T& w) {
	if (w.unit())
	    return *this;
	typename factorlist<T>::iterator i;
	for(i = begin(); i!= end() && *i<w; ++i);
	if (i!= end() && *i==w)
	    *i *= w;
	else
	    insert(i,w);
	return *this;
    }
    factorlist<T>& operator=(const T& w) {
	clear();
	if (!w.unit())
	    push_front(w);
	return *this;
    }

/*      bool operator==(const factorlist<T>& W) const { */
/*  	factorlist<T>::const_iterator i = begin(); */
/*  	factorlist<T>::const_iterator wi = W.begin(); */
/*  	for (; i != end() && wi != W.end(); ++i, ++wi) */
/*  	    if (!(*i == *wi && i->exp() == wi->exp())) */
/*  		return false; */
/*  	return i == end() && wi == W.end(); */
/*      } */

    factorlist<T>& operator*=(factorlist<T>& W) { 
	typename factorlist<T>::iterator i = begin();
	typename factorlist<T>::iterator wi = W.begin();
	while (!W.empty()) {
	    for(; i != end() && *i < *wi; ++i);
	    if (i == end()) {
		splice(i,W);
		break;
	    }
	    for(; wi != W.end() && *wi < *i; ++wi);
	    if (wi != W.begin()) {
		splice(i,W,W.begin(),wi);
		if (W.empty())
		    break;
		wi = W.begin();
	    }
	    if (*wi == *i) {
		*i *= *wi;
		W.erase(wi);
		wi = W.begin();
		if (i->unit())
		    i = erase(i);
		else
		    ++i;
	    }
	    else
		++i;
	}
	return *this;
    }

    factorlist<T>& operator^=(int q) {
	assert(q>=0);
	if (q == 0) {
	    clear();
	    return *this;
	}
	for(typename factorlist<T>::iterator i = begin(); i != end(); i++) {
	    *i ^= q;
	}
	return *this;
    }

    void move(factorlist<T>& W) {
	clear();
	splice(end(),W);
    }
};

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

for expressions of the form "b + mt", used for the exponents of the
  z terms

for more general exponents it might be possible to parameterize all 
  uses of lin

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

class lin {
private:
    int b;
    int m;
public:
    lin(int b1, int m1) {b = b1; m = m1;}
    lin() {b = 0; m = 0;}
    int tpart() const {return m;}
    int ipart() const {return b;}

// evaluation as a poly in t
    poly& tpoly() const;

    lin& operator=(const int b1) {
	b = b1;
	m = 0;
	return *this;
    }
    lin& operator+=(const lin& z) {
	b += z.b;
	m += z.m;
	return *this;
    }
    lin& operator+=(const int z) {
	b += z;
	return *this;
    }
    lin& operator-=(const lin& z) {
	b -= z.b;
	m -= z.m;
	return *this;
    }
    lin& operator-=(const int z) {
	b -= z;
	return *this;
    }
    lin& operator*= (int w) {
	b *= w;
	m *= w;
	return *this;
    }

// this is lexicographic order -- in effect, treat "t" as infinite and positive

    bool operator< (const lin& w) const {
	if (m != w.m)
	    return m < w.m;
	return b < w.b;
    }
    bool operator< (int w) const {
	return m < 0 || (m == 0 && b < w);
    }
    bool operator== (const lin& w) const {
	return m == w.m && b == w.b;
    }
    bool operator== (const int y) const {
	return m == 0 && b == y;
    }

    friend ostream& operator<<(ostream& f, const lin& z);
};

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

the binomial coefficient [n : k]

these arise form integration, so n may involve the exponents of the z's,
  so n must be of type lin

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

class binomial {
private:
    lin n;
    int k;
public:
    binomial(const lin& n1, int k1) :n(n1) {k = k1;}
    binomial(int b1, int m1, int k1) :n(b1,m1) {k = k1;}
    binomial() {k = 0;}
    const lin& top() const {return n;}
    int bottom() const {return k;}

// evaluation as a poly in t
    poly& tpoly() const;

// the topcoefficient of tpoly -- returns false and does nothing if
// this is a constant in t
    bool topcoeff(RAT& v) const;

// bombs if this depends on t
    INT& eval(INT& v) const;

// the ordering is used only by factorlist

    bool operator< (const binomial& w) const {
	if (n < w.n)
	    return true;
	if (w.n < n)
	    return false;
	return k < w.k;
    }
    bool operator== (const binomial& w) const {
	return n == w.n && k == w.k;
    }

    friend ostream& operator<<(ostream& f, const binomial& c);
};

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

a bfactor is just a binomial to a power

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

class bfactor {
private:
    binomial bc;
    int p;
public:
    bfactor(lin& n1, int k1, int p1) :bc(n1,k1) {p = p1; assert(p);}
    bfactor(int n1, int k1, int p1) :bc(n1,0,k1) {p = p1; assert(p);}
    bfactor(int b1, int m1, int k1, int p1) :bc(b1,m1,k1) {p = p1; assert(p);}
    bfactor() :bc(0,0,0) {p=1;}
    const binomial& base() const {return bc;}
    int exp() const {return p;}

// used for a hash key by upoly

    uint cache_key() const {
	assert(-127 < base().top().ipart() &&
	       base().top().ipart() <= 127);
	return (((unsigned)base().top().ipart())<<24) |
	    (base().top().tpart()<<16) |
	    (base().bottom() << 8) |
	    (exp());
    }

private:
    const poly& make_tpoly() const;
public:
    const poly& tpoly() const;

// bombs if this depends on t
    INT& eval(INT& v) const {
	return bc.eval(v) ^= p;
    }
    bool topcoeff(RAT& v) const {
	if (!bc.topcoeff(v))
	    return false;
	v ^= p;
	return true;
    }

// these are for use by factorlist
    bool operator< (const bfactor& w) const {
	return bc < w.bc;
    }
    bool operator== (const bfactor& w) const {
	return bc == w.bc;
    }
/*      bool unit() const {return p == 0;} */
    bool unit() const {return p == 0 || bc.top() == bc.bottom();}
    bfactor& operator*=(const bfactor& w) {
	p += w.p;
	return *this;
    }
    bfactor& operator^= (int q) {
	p *= q;
	return *this;
    }
	    
    friend ostream& operator<<(ostream& f, const bfactor& c);

// only used to set up hash keys
    string varstr() const {
	ostringstream os;
	os << bc.top() << ";" << bc.bottom() << ";" << p;
	return os.str();
    }
};

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

a zfactor represents a factor in the numerator, of the form z_k^p

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

class zfactor {
private:
    int v;
    lin p;
public:
    zfactor () {}
    zfactor(int v1, int p1) :p(p1,0) {
	v = v1;
    }
    zfactor(int v1, int b1, int m1) :p(b1,m1) {
	v = v1;
    }
    zfactor(int v1, const lin& p1) :p(p1) {
	v = v1;
    }
    int var() const {return v;}
    const lin& exp() const {return p;}

// these are for use by factorlist
    bool operator< (const zfactor& w) const {
	 return v < w.v;
    }
    bool operator== (const zfactor& w) const {
	return v == w.v;
    }
// we don't erase (z_k)^0 factors since integrate needs to see them
    bool unit() const {return false;}
    zfactor& operator*=(const zfactor& w) {
	p += w.p;
	return *this;
    }
    zfactor& operator^= (int q) {
	p *= q;
	return *this;
    }

/********************
  maxzdiff is used by term::residue and also by statterm::evaluate

  it calculates the maximum number of times that this zfactor should be 
  differentiated (with an a priori upper bound of q)

  this is the only place where TopCoefficient is used --
  if TopCoefficient is true then maxzdiff is 0 if this zfactor has an
  exponent not containing t.
  in this case, differentiating the zfactor will lower the t degree of
  the result, so we will not be able to produce top degree in t
********************/

    int maxzdiff(int q, bool TopCoefficient) const {
	int maxq;
	if (exp().tpart())
	    maxq = q;
	else if (TopCoefficient)
	    maxq = 0;
	else {
	    maxq = exp().ipart();
	    if (maxq < 0 || maxq > q)
		maxq = q;
	}
	return maxq;
    }

    friend ostream& operator<<(ostream& f, const zfactor& c);
};    

typedef factorlist<bfactor>::iterator BI;
typedef factorlist<bfactor>::const_iterator CBI;
typedef factorlist<zfactor>::iterator ZI;
typedef factorlist<zfactor>::const_iterator CZI;
typedef factorlist<zfactor>::reverse_iterator RZI;
typedef factorlist<zfactor>::const_reverse_iterator CRZI;

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

a zzfactor represents a factor in the denominator, of the form 
(z_k - z_j)^p.  p is stored as a negative integer.

most routines assume that k>j in a zzfactor; normalorder() is used to 
ensure this

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

class zzfactor {
private:
    int var1,var2;
    int p;
public:
    zzfactor () {}
    zzfactor(int v1, int v2, int p1) {
	var1 = v1;
	var2 = v2;
	p = p1;
    }
    int var(int which) const {return which == 1 ? var1: var2;}
    int exp() const {return p;}

// interchange vars if necessary; return the change in sign
    int normalorder() {
	if (var1 > var2)
	    return 1;
	int var = var1;
	var1 = var2;
	var2 = var;
	if(ODD(p))
	    return -1;
	else
	    return 1;
    }

// these are for use by factorlist
    bool operator< (const zzfactor& w) const {
	if (var1 < w.var1)
	    return true;
	if (w.var1 < var1)
	    return false;
	return var2 < w.var2;
    }
    bool operator== (const zzfactor& w) const {
	return var1 == w.var1 && var2 == w.var2;
    }
    bool unit() const {return p == 0;}
    zzfactor& operator*=(const zzfactor& w) {
	p += w.p;
	return *this;
    }
    zzfactor& operator^= (int q) {
	assert(q>0);
	p *= q;
	return *this;
    }

    friend ostream& operator<<(ostream& f, const zzfactor& c);
};    

typedef factorlist<zzfactor>::iterator ZZI;
typedef factorlist<zzfactor>::const_iterator CZZI;
typedef factorlist<bfactor> bprod;
typedef factorlist<zfactor> zprod;
typedef factorlist<zzfactor> zzprod;

enum METHOD {M_EXT, M_INT, M_BRANCHCOUNT, M_QUERY, M_TOPCOEFF};

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

a term represents an integrand -- a product of a scalar, bfactors, zfactors and
zzfactors.

terms keep each zfactor (possibly with 0 exponent) until integration on that
variable.  also, all variables that appear in the zzfactors appear as zfactors

the static initialize member sets up some static variables in algebra.cpp
********************/

class term {
    friend class statterm;
    friend class tterm;
private:
    int s;
    bprod B;
    zprod Z;
    zzprod ZZ;
public:
    static CTR numresidue;
    static CTR numint;
    static CTR numzeroint;
    static CTR numsimpleint;
    static CTR numdoubleint;
    static CTR nummult;
public:
    static void initialize(int n, int topdeg, bool tracing, 
			   bool factoring, bool zerocheck, METHOD method);
public:
    term() {s=0;}
    int scomp() const {return s;}
    int numz() const {return Z.size();}

// misleading name: constant() really means "completely integrated"
    int constant() const {return Z.empty();}
    int zzdeg(int v) const {
	int deg = 0;
	for(CZZI izz = ZZ.begin(); izz != ZZ.end(); ++izz)
	    if (v == izz->var(1) || v == izz->var(2))
		deg += izz->exp();
	return deg;
    }
    int zzdeg(int v1, int v2) const {
	if (v1<v2)
	    return zzdeg(v2,v1);
	for(CZZI izz = ZZ.begin(); izz != ZZ.end(); ++izz)
	    if (v1 == izz->var(1) && v2 == izz->var(2))
		return izz->exp();
	return 0;
    }
    term& operator=(int s1) {
	s=s1;
	B.clear();
	Z.clear();
	ZZ.clear();
	return *this;
    }
    term& operator=(const bfactor& f) {
	Z.clear();
	ZZ.clear();
	s=1;
	B = f;
	return *this;
    }
    term& operator=(const zfactor& f) {
	B.clear();
	ZZ.clear();
	s=1;
	Z = f;
	return *this;
    }
    term& operator=(const zzfactor& f) {
	B.clear();
	Z.clear();
	s=1;
	ZZ = f;
	return *this;
    }
    term& operator*= (int w) {
	if(w)
	    s *= w;
	else
	    *this = 0;
	return *this;
    }
    term& operator*= (term& w) {
	if (w.s) {
	    s *= w.s;
	    B *= w.B;
	    Z *= w.Z;
	    ZZ *= w.ZZ;
	}
	else
	    *this = 0;
	return *this;
    }
    term operator*= (zfactor& w) {
	Z *= w;
	return *this;
    }
    term operator*= (zzfactor& w) {
	ZZ *= w;
	return *this;
    }
    term& operator^= (int q) {
	assert(q>=0);
	int t = s;
	s = 1;
	for (int i=0; i<q; i++)
	    s *= t;
	B ^= q;
	Z ^= q;
	ZZ ^= q;
	return *this;
    }
    void move(term& w) {
	B.move(w.B);
	Z.move(w.Z);
	ZZ.move(w.ZZ);
	s = w.s;
    }

    void topcoeff(RAT& v) {
	v = s;
	RAT x;
	INT y;
	for(CBI ib = B.begin(); ib != B.end(); ++ib)
	    if (ib->topcoeff(x))
		v *= x;
	    else
		v *= ib->eval(y);
    }
    INT& evalcoeff(INT& v) const {
	INT f;
	v = s;
	for(CBI ib = B.begin(); ib != B.end(); ++ib)
	    v *= ib->eval(f);
	return v;
    }
    int topdeg() const {
	int d = 0;
	for(CBI ib = B.begin(); ib != B.end(); ++ib)
	    if(ib->base().top().tpart())
		d += ib->exp()*ib->base().bottom();
	return d;
    }
    poly& tpoly() const;

    generator* integrate1();

private:
    void splitoff(int v, term & factor);
    generator* zero_integral();
    generator* simple_integral();
    generator* double_integral();
    generator* ext_integral(CZI startz, int multiplicity);
    generator* int_integral(CZI startz, int multiplicity);
    bool zero_res_at_zero (ZI iz) const;
    bool zero_res_at_infinity (ZI iz) const;
    expr& residue(int r, int var, int pole, expr& result);
    expr& zzresidue(int q, int var, int pole, expr& result);
    term& zresidue1(int r, int var, int pole);
    term& zzresidue1(int r, int var, int pole);

public:
    friend ostream& operator<<(ostream& f, const term& T);
};

/********************
  a tterm is a term with the constant bfactors converted to INTs
  it is a friend of the term class; this is used only in the conversion
    function operator= which iterates over the B member of a term
  this can be avoided by providing a term member function that provides
    B as a const reference
********************/

class tterm {
private:
    bprod B;
    INT S;
public:
    tterm() {}
    ~tterm() {}
    const INT& Scomp() const {return S;}
    tterm& operator= (const term& T) {
	B.clear();
	S = T.scomp();
	INT V;
	for(CBI ib = T.B.begin(); ib != T.B.end(); ++ib)
	    if (ib->base().top().tpart())
		B *= *ib;
	    else
		S *= ib->eval(V);
	return *this;
    }
    void topcoeff(RAT& v) const {
	v = S;
	RAT x;
	for(CBI ib = B.begin(); ib != B.end(); ++ib) {
	    ib->topcoeff(x);
	    v *= x;
	}
    }
    int topdeg() const {
	int d = 0;
	for(CBI ib = B.begin(); ib != B.end(); ++ib) {
	    d += ib->exp()*ib->base().bottom();
	}
	return d;
    }
    bool zero() const {return sgn(S) == 0;}
    tterm& operator+= (const tterm& T) {
	if (zero())
	    *this = T;
	else {
	    assert(B == T.B);
	    S += T.S;
	    if (zero())
		B.clear();
	}
	return *this;
    }
    poly& tpoly() const;

    friend ostream& operator<<(ostream& f, const tterm& T);
    string varstr() const {
	ostringstream os;
	if (B.empty())
	    os << 1;
	else
	    for (CBI ib = B.begin(); ib != B.end(); ++ib)
		os << *ib;
	return os.str();
    }
};

typedef list<term>::iterator TI;
typedef list<term>::const_iterator CTI;

typedef list<term> tsum;

/********************
  an expr is just a sum of terms
********************/

class expr {
private:
    tsum T;
public:
    expr() {}
    bool empty() const {return T.empty();}
    int numterms() const {return T.size();}
    void getfirst(term& t) {
	assert(!T.empty());
	t.move(T.front());
	T.pop_front();
    }
    void addfirst(expr& W) {
	T.splice(T.begin(),W.T);
    }
    expr& operator= (const term &w) {
	T.clear();
	if (w.scomp() == 0)
	     return *this;
	term t;
	T.push_back(t);
	T.back() = w;
	return *this;
    }
    expr& operator= (int w) {
	T.clear();
	if (w == 0)
	    return *this;
	term t;
	T.push_back(t);
	T.back() = w;
	
	return *this;
    }
    expr& operator+= (term &w) {
	if (w.scomp()) {
	    term t;
	    T.push_back(t);
	    T.back().move(w);
	}
	return *this;
    }
    expr& operator+= (expr& w) {
	T.splice(T.end(),w.T);
	return *this;
    }

    expr& operator*= (term& w) {
	if (empty())
	    return *this;
	if (w.scomp() == 0) {
	    *this = 0;
	    return *this;
	}
	term t;
	TI j = T.end();
	--j;
	for (TI i=T.begin(); i!=j; ++i) {
	    t = w;
	    *i *= t;
	}
	*j *= w;
	return *this;
    }
    expr& operator*= (expr& w);

// this is not currently used, and should probably be rewritten 
// as an operator^=, if feasible
/*      expr& pow(expr& E, int p); */
/*  private: */
/*      expr& makepow(expr& E, int m, term& t, TI it); */
    friend ostream& operator<<(ostream& f, const expr& E);
};

#endif