using namespace std;
#include "bignum.h"
#include "count.h"
#include "calc.h"
#include "poly.h"
#include "algebra.h"
#include "gen.h"
#include "upoly.h"
#include "statterm.h"
#include "utilio.h"

upoly U;

ctr upoly::numsetcache;
CTR upoly::numgetcache;

statterm S;
expr scratchE;

bool Tracing = false;
bool Factoring = false;
bool TopCoefficient = false;
bool ZeroCheck = false;
METHOD Method = M_BRANCHCOUNT;

/********************
  IO operators
********************/

ostream& put_sym(ostream& f, const int m, char * sym, bool showsign) {
    char* sign = "";
    if (m==0)
	return showsign ? f : f << 0;
    if (m>0) {
	if (showsign)
	    sign = "+";
	if (m==1)
	    return f << sign << sym;
	else
	    return f << sign << m << sym;
    }
    if (m==-1)
	return f << "-" << sym;
    return f << m << sym;
}

ostream& operator<<(ostream& f, const lin& z) {
    if (z.b == 0)
	return put_sym(f,z.tpart(),"t",false);
    return put_sym(f<<z.b,z.tpart(),"t",true);
}

ostream& operator<<(ostream& f, const binomial& c) {
    f << "[" << c.top() << ":" << c.bottom() << "]";
    return f;
}

ostream& operator<<(ostream& f, const bfactor& c) {
    f << c.base();
    if (c.exp()!=1)
	if (c.exp() >= 0)
	    f << "^" << c.exp();
	else
	    f << "^{" << c.exp() << "}";
    return f;
}

ostream& operator<<(ostream& f, const zfactor& c) {
    char * lpar = "(";
    char * rpar = ")";

    f << lpar << "z" << c.var() << rpar;
    if (!(c.exp()==1))
	if (c.exp().ipart() > 0 && c.exp().tpart() == 0)
	    f << "^" << c.exp();
	else
	    f << "^{" << c.exp() << "}";
    return f;
}

ostream& operator<<(ostream& f, const zzfactor& c) {
    char * lpar = "(";
    char * rpar = ")";
    if (c.var(1))
	if (c.var(2))
	    f << lpar << "z" << c.var(1) << "-z" << c.var(2) << rpar;
	else
	    f << lpar << "z" << c.var(1) << rpar;
    else
	if (c.var(2))
	    f << lpar << "-z" << c.var(2) << rpar;
	else
	    f << 0;
    if (c.exp()!=1)
	if (c.exp() >= 0)
	    f << "^" << c.exp();
	else
	    f << "^{" << c.exp() << "}";
    return f;
}

ostream& operator<<(ostream& f, const term& T) {
    if (T.s == 0)
	return f << 0;
    if (T.B.empty() && T.Z.empty() && T.ZZ.empty()) {
	return f << T.s;
    }
    if (T.s == -1)
	f << "-";
    else if (T.s != 1)
	f << T.s;
    for (CBI ib = T.B.begin(); ib != T.B.end(); ++ib)
	f << *ib;
    for (CZI iz = T.Z.begin(); iz != T.Z.end(); ++iz)
	f << *iz;
    for (CZZI izz = T.ZZ.begin(); izz != T.ZZ.end(); ++izz)
	f << *izz;
    return f;
}

ostream& operator<<(ostream& f, const tterm& T) {
    if (sgn(T.S) == 0)
	return f << 0;
    if (T.B.empty())
	return f << T.S;
    if (T.S == -1)
	f << "-";
    else if (T.S != 1)
	f << T.S;
    for (CBI ib = T.B.begin(); ib != T.B.end(); ++ib)
	f << *ib;
    return f;
}

ostream& operator<<(ostream& f, const expr& E) {
    if (E.T.empty())
	return f << "0\n";
    int termnum = 0;
    for(CTI i = E.T.begin(); i != E.T.end(); ++i) {
	if (termnum++)
	    if (i->scomp() >= 0)
		f << "+";
	f << *i << "\n";
    }
    return f;
}

/********************
  tpoly calculations
********************/

poly& lin::tpoly() const {
    poly& Q = U.scratch(0);
    Q.set_numer(0,ipart());
    Q.set_numer(1,tpart());
    return Q;
}

poly& binomial::tpoly() const {
    poly& F = top().tpoly();
    poly& Q = U.scratch(1);
    Q = F;
    for(int j = 2; j <= bottom(); j++) {
	F.set_numer(0,F.get_numer(0)-1);
	Q *= F;
    }
    Q /= C.Factorial(bottom());
    return Q;
}

const poly& bfactor::make_tpoly() const {
    poly& P = base().tpoly();
    P ^= exp();
    return P;
}

const poly& bfactor::tpoly() const {
    if (base().bottom()*exp() <= 1) 
	return make_tpoly();
    else {
	uint key = cache_key();
	const poly& P = U.get_cache_poly(key);
	if (P.deg() > 0)
	    return P;
	else
	    return U.set_cache_poly(key,make_tpoly());
    }
}

poly& term::tpoly() const {
    poly& P = U.scratch(3);
    P = 1;
    INT v = s;
    INT y;
    for(CBI ib = B.begin(); ib != B.end(); ++ib)
	if(ib->base().top().tpart())
	    P *= ib->tpoly();
	else
	    v *= ib->eval(y);

    P *= v;
    return P;
}

poly& tterm::tpoly() const {
    poly& P = U.scratch(3);
    P = 1;
    INT v = S;
    INT y;
    for(CBI ib = B.begin(); ib != B.end(); ++ib)
	if(ib->base().top().tpart())
	    P *= ib->tpoly();
	else
	    v *= ib->eval(y);
	
    P *= v;
    return P;
}

/********************
  mult and ops on exprs
********************/

expr& expr::operator*= (expr& w) {
    if (w.empty()) {
	*this = 0;
	return *this;
    }
    term t;
    w.getfirst(t);
    if (w.empty()) {
	*this *= t;
	return *this;
    }
    expr F;
    expr E0;
    E0 = *this;
    *this *= t;
    F = E0;
    while (true) {
	w.getfirst(t);
	if (!w.empty()) {
	    F = E0;
	    F *= t;
	    *this += F;
	}
	else {
	    E0 *= t;
	    *this += E0;
	    break;
	}
    }
    return *this;
}

//  expr& expr::makepow(expr& E, int m, term& t, TI it) {
//      if (m == 0) {
//  	E += t;
//  	return E;
//      }
//      term tc;
//      if (it == T.begin()) {
//  	tc = *it;
//  	tc ^= m;
//  	t *= tc;
//  	E += t;
//  	return E;
//      }
//      TI nit = it;
//      --nit;
//      int c = 1;
//      term t1;
//      t1 = t;
//      makepow(E,m,t1,nit);
//      for (int q = 1; q <= m; q++) {
//  	tc = *it;
//  	tc ^= q;
//  	t1 = t;
//  	t1 *= tc;
//  	c *= m - q + 1;
//  	c /= q;
//  	t1 *= c;
//  	makepow(E,m-q,t1,nit);
//      }
//      return E;
//  }

//  expr& expr::pow(expr& E, int p) {
//      assert(p>=0);
//      if (p == 0) {
//  	E = 1;
//  	return E;
//      }
//      E = 0;
//      if (T.empty())
//  	return E;
//      term t;
//      t = 1;
//      TI it = T.end();
//      --it;
//      return makepow(E,p,t,it);
//  }

/********************
  coefficient calculations
********************/

INT& binomial::eval(INT& v) const {
    assert(top().tpart() == 0);
    return C.Ibin(v,top().ipart(),bottom());
}

bool binomial::topcoeff(RAT& v) const {
    if (top().tpart() == 0)
	return false;
    INT V;
    V = top().tpart();
    if (top().tpart() != 1)
	V ^= bottom();
    fraction(v,V,C.Factorial(bottom()));
    return true;
}
 
/********************
  integration
********************/

// splits off the factors in *this that involve var
void term::splitoff(int v, term& factor) {
    factor = 1;
    for (ZI iz = Z.begin(); iz != Z.end(); ++iz) 
        if (iz->var() == v) {
            factor.Z.splice(factor.Z.end(),Z,iz);
            break;
        }
    ZZI jzz;
    for (ZZI izz = ZZ.begin(); izz != ZZ.end(); ) 
	if (izz->var(1) == v || izz->var(2) == v) {
	    jzz = izz;
	    ++izz;
	    factor.ZZ.splice(factor.ZZ.end(),ZZ,jzz);
	}
	else
	    ++izz;
    if (Tracing) cerr << "after split on z" << v << ": \n" << 
		       "    " << *this << "\n" <<
		       "    " << factor << "\n";
}

/********************
  the various residue functions operate on the term t = *this with exponent 
  r, variable z and pole b.  in fact, the function calculates the residue of 
  t/(z-b)^(r+1) at b, so this is calculated as the rth derivative of t, 
  evaluated at b, divided by r!.

  using this definition, the residue of a product is the
  Leibnitz sum of the products of the residues, and this is used in
********************/

// only called if exactly one zfactor, which must contain var (and not pole),
// no zzfactors or bfactors,
// so that the rth derivative is not 0
// replaces *this with its residue

term& term::zresidue1(int r, int var, int pole) {
    lin p = Z.front().exp();
    if (r>0)
	if (p < 0) {
	    B = bfactor(-p.ipart()+r-1,-p.tpart(),r,1);
	    if (ODD(r))
		s = -s;
	}
	else
	    B = bfactor(p,r,1);
    p -= r;
    Z.front() = zfactor(pole,p);
    return *this;
}

// only called if exactly one zzfactor, which must contain var (and not pole)
// no zfactor or bfactor; assumes that the zzfactor has negative exponent
// replaces *this with its residue

term& term::zzresidue1(int r, int var, int pole) {
    int p = ZZ.front().exp();
    if (r>0) {
	B = bfactor(-p+r-1,r,1);
	if (ODD(r))
	    s = -s;
    }
    p -= r;
    if (ZZ.front().var(1) == var)
	ZZ.front() = zzfactor(pole,ZZ.front().var(2),p);
    else {
	if (ODD(r))
	    s *= -1; // chain rule
	ZZ.front() = zzfactor(ZZ.front().var(1),pole,p);
    }
    s *= ZZ.front().normalorder();
    return *this;
}

// assumes no Z factors and at least one ZZ factor; 
// all ZZ factors contain var (but not pole)
// destroys *this, provides answer in result

expr& term::zzresidue(int q, int var, int pole, expr& result) {
    numresidue++;
    if (Tracing) cerr << "finding residue(" << numresidue << 
			    ") wrt z" << var << " at z" << 
			    pole << " of\n" << "    " << *this << "\n";
    term f;
    f = 1;
    f.ZZ.splice(f.ZZ.end(),ZZ,ZZ.begin());

    if (ZZ.empty()) {
// Get any existing binomials and scalars in *this
	f *= *this;
	return result = f.zzresidue1(q,var,pole);
    }

    term tc;
    expr E;
    for(int r=0; r<=q; r++) {
	tc = *this;
	tc.zzresidue(q-r,var,pole,E);
	tc = f;
	E *= tc.zzresidue1(r,var,pole);
	result += E;
    }
    return result;
}

// assumes there is exactly one zfactor
// destroys *this, provides answer in result

expr& term::residue(int q, int var, int pole, expr& result) {
    
    assert(!Z.empty());
    numresidue++;
    
    if (Tracing) cerr << "finding residue(" << numresidue << 
			    ") wrt z" << var << " at z" << 
			    pole << " of\n" << "    " << *this << "\n";

    term f;
    f = 1;
    f.Z.splice(f.Z.end(),Z,Z.begin());

    int maxr = f.Z.front().maxzdiff(q,TopCoefficient);
    
    if (ZZ.empty()) {
	if (q > maxr)
	    result = 0;
	else {
// Get any existing binomials and scalars in *this
	    f *= *this;
	    f.zresidue1(q,var,pole);
	    result = f;
	}
	return result;
    }

    term tc;
    expr E;
    for(int r=0; r<=maxr; r++) {
	tc = *this;
	tc.zzresidue(q-r,var,pole,E);
	tc = f;
	E *= tc.zresidue1(r,var,pole);
	result += E;
    }
    return result;
}

inline bool term::zero_res_at_zero (ZI iz) const {
    return !(iz->exp() < 0);
}

inline bool term::zero_res_at_infinity (ZI iz) const {
    return iz->exp().tpart() < 0 ||
	iz->exp().tpart() == 0 && iz->exp().ipart() + zzdeg(iz->var()) < -1;
}

// choose one variable and direction and integrate
generator* term::integrate1() {
    numint++;
    if (Tracing) cerr << "integrating (" << numint << "): " << 
			   *this << "\n";
    switch(numz()) {
    case 0:
	if (Tracing) cerr << "constant term; integral is 0\n";
	return zero_integral();
    case 1:
	return simple_integral();
    case 2:
	return double_integral();
    }

    ZI ez,iz;

    switch (Method) {
    case M_EXT:
	ez = Z.end();
	--ez;
	if (zero_res_at_zero(ez)) {
	    if (Tracing) cerr << "analytic in z" << ez->var() << 
				   " inside circle r" << ez->var() << "\n";
	    return zero_integral();
	}
	assert(zero_res_at_infinity(ez));
	return ext_integral(ez,1);
    case M_INT:
	iz = Z.begin();
	if (zero_res_at_infinity(iz)) {
	    if (Tracing) cerr << "order < -1 in z" << iz->var() << 
				   " outside circle r" << iz->var() << "\n";
	    return zero_integral();
	}
	assert(zero_res_at_zero(iz));
	return int_integral(iz,1);
    case M_QUERY:
	ez = Z.end();
	--ez;
	if (zero_res_at_zero(ez)) {
	    if (Tracing) cerr << "analytic in z" << ez->var() << 
				   " inside circle r" << ez->var() << "\n";
	    return zero_integral();
	}
	iz = Z.begin();
	if (zero_res_at_infinity(iz)) {
	    if (Tracing) cerr << "order < -1 in z" << iz->var() << 
				   " outside circle r" << iz->var() << "\n";
	    return zero_integral();
	}
	{
	    int qvar = 0;
	    char qdir;
	    int qmult;
	    char errch;
	    CZI qz;
	    while (qvar == 0) {
		cout << "Integrating " << *this << "\n";
		S.evaluate(*this);
		cout << "    Evaluations:\n";
		S.dumpeval(cout);
		if (ZeroCheck && S.zeroz) {
		    cout << "(integral is zero) ";
		    if(skiptoeol(cin) == EOF) {
			cout << "\n";
			exit(1);
		    }
		    return zero_integral();
		}
		cout << "Enter varnum dir [mult]: ";

		if (!parseline(cin,qvar,qdir,qmult,errch)) {
		    cout << "\n";
		    exit(1);
		}
		if(errch) {
		    cout << "can\'t interpret \'" << errch
			 << "\' (" << (int)errch << ")\n";
		    qvar = 0;
		    continue;
		}
		if(!qmult)
		    qmult = 1;
		if(!qvar || !qdir) {
		    cout << "must specify varnum and dir\n";
		    qvar = 0;
		    continue;
		}
		if(!S.ok_int(qvar,qdir,qmult,qz)) {
		    cout << "illegal choice\n";
		    qvar = 0;
		    continue;
		}
		if (qdir == 'e') 
		    return ext_integral(qz,qmult);
		else // qdir == 'i'
		    return int_integral(qz,qmult);
	    }
	}
	break;
    case M_TOPCOEFF:
    case M_BRANCHCOUNT:
	ez = Z.end();
	--ez;
	if (zero_res_at_zero(ez)) {
	    if (Tracing) cerr << "analytic in z" << ez->var() << 
				   " inside circle r" << ez->var() << "\n";
	    return zero_integral();
	}
	iz = Z.begin();
	if (zero_res_at_infinity(iz)) {
	    if (Tracing) cerr << "order < -1 in z" << iz->var() << 
				   " outside circle r" << iz->var() << "\n";
	    return zero_integral();
	}
	S.evaluate(*this);
	
	if (Tracing) S.dumpeval(cerr);

	if(ZeroCheck && S.zeroz)
	    return zero_integral();
	if(S.bestdir == 'e') 
	    return ext_integral(S.bestz, S.bestmult);
	else
	    return int_integral(S.bestz, S.bestmult);

    default:
	assert(0);
    }
    return 0; // not reached
}

// zero integral, for various reasons
// a separate function just to centralize treatment
inline generator* term::zero_integral() {
    numzeroint++;
    return 0;
}    

// only one variable, so no zz terms
generator* term::simple_integral() {
    if (Tracing) cerr << "simple integral of " << *this << ":\n";
    numsimpleint++;
    if (Z.begin()->exp() == -1) {
	Z.clear();
	scratchE = *this;
	if (Tracing) cerr << "    " << scratchE;
	expgen* GE = new expgen;
	GE->initialize(scratchE);
	return GE;
    }
    else
	return zero_integral();
}

// two variables, so at most one zz term
generator* term::double_integral() {
    if (Tracing) cerr << "double integral of " << *this << ":\n";
    numdoubleint++;
    ZI iz1 = Z.begin();
    ZI iz2 = iz1;
    ++iz2;
    if (ZZ.empty()) {
	if (!(iz1->exp() == -1 && iz1->exp() == -1))
	    return zero_integral();
	Z.erase(iz1);
	return simple_integral();
    }
    if (iz2->exp() < 0) {
	int q = ZZ.begin()->exp();
	B *= bfactor(iz1->exp().ipart(),iz1->exp().tpart(),-q-1,1);
	if (ODD(q))
	    s *= -1;
	ZZ.clear();
	Z.clear();
	scratchE = *this;
	if (Tracing) cerr << "    " << scratchE;
	expgen* GE = new expgen;
	GE->initialize(scratchE);
	return GE;
    }
    else
	return zero_integral();
}

// the residue at infinity is 0 and the residue at zero is not 0
// multiplicity entries starting at iz can be combined
generator* term::ext_integral(CZI startz, int multiplicity) {
    term f;
// initializing to silence complaints about var possibly being used
// unitialized (can't happen, of course, since multiplicity >= 1)
    int var = 0; 

    if(multiplicity > 1)
	nummult++;
    for (int m=1; m<=multiplicity; m++) {
	var = startz->var();
	++startz;
	splitoff(var,f);
    }
	
// I don't think this can happen ...
    if (f.ZZ.empty()) {
	ZI iz = f.Z.begin();
	if (Tracing) cerr << "just integrating " <<
				*iz << "; result is ";
	if (iz->exp() == -1) {
	    if (Tracing) cerr << "1\n";
	    scratchE = *this;
	    expgen* GE = new expgen;
	    GE->initialize(scratchE);
	    return GE;
	}
	else {
	    if (Tracing) cerr << "0\n";
	    return zero_integral();
	}
    }
    
    scratchE = 0;
    int pole;
    expr result;
    int q;
    term g;

// summing over exterior poles:
    f.s *= -1;

    for (ZZI i = f.ZZ.begin(); i != f.ZZ.end(); ++i) {
// assuming normal order, this means var(1) > var, so not an exterior pole
	if (i->var(2) == var)
	    continue;
	pole = i->var(2);
	q = -i->exp();

	g = f;

	ZZI j = g.ZZ.begin();
	for (; j->var(2) != pole; ++j);
	g.ZZ.erase(j);

	g.residue(q-1,var,pole,result);
	if (Tracing) cerr << "result:\n" << result;
	scratchE += result;
    }
    if (multiplicity == 1) {
	scratchE *= *this;
	if (Tracing) cerr << "adding to integral:\n" << scratchE;
	expgen* GE = new expgen;
	GE->initialize(scratchE);
	return GE;
    }
    else {
	if (Tracing) cerr << "adding to integral:\n" << *this << "\n* (\n" 
			    << scratchE << ")^" << multiplicity << " \n";
	powgen* GP = new powgen;
	GP->initialize(*this,scratchE,multiplicity);
	return GP;
    }
}

// the residue at infinity is 0 and the residue at infinity is not 0
// multiplicity entries starting at iz can be combined
generator* term::int_integral(CZI startz, int multiplicity) {
    term f;
// initializing to silence complaints about var possibly being used
// unitialized (can't happen, of course, since multiplicity >= 1)
    int var = 0; 

    if(multiplicity > 1)
	nummult++;
    for (int m=1; m<=multiplicity; m++) {
	var = startz->var();
	++startz;
	splitoff(var,f);
    }
	
// I don't think this can happen ...
    if (f.ZZ.empty()) {
	ZI iz = f.Z.begin();
	if (Tracing) cerr << "just integrating " <<
				*iz << "; result is ";
	if (iz->exp() == -1) {
	    if (Tracing) cerr << "1\n";
	    scratchE = *this;
	    expgen* GE = new expgen;
	    GE->initialize(scratchE);
	    return GE;
	}
	else {
	    if (Tracing) cerr << "0\n";
	    return zero_integral();
	}
    }
    
    scratchE = 0;
    int pole;
    expr result;
    int q;
    term g;

    for (ZZI i = f.ZZ.begin(); i != f.ZZ.end(); ++i) {
// assuming normal order, this means var(2) < var, so not an exterior pole
	if (i->var(1) == var)
	    continue;
	pole = i->var(1);
	q = -i->exp();

	g = f;
// since zzfactor is (pole - var)^q, not (var - pole)^q
  	if (ODD(q))
  	    g *= -1;

	ZZI j = g.ZZ.begin();
	for (; j->var(1) != pole; ++j);
	g.ZZ.erase(j);

	g.residue(q-1,var,pole,result);
	if (Tracing) cerr << "result:\n" << result;
	scratchE += result;
    }
    if (multiplicity == 1) {
	scratchE *= *this;
	if (Tracing) cerr << "adding to integral:\n" << scratchE;
	expgen* GE = new expgen;
	GE->initialize(scratchE);
	return GE;
    }
    else {
	if (Tracing) cerr << "adding to integral:\n" << *this << "\n* (\n" 
			    << scratchE << ")^" << multiplicity << " \n";
	powgen* GP = new powgen;
	GP->initialize(*this,scratchE,multiplicity);
	return GP;
    }
}

CTR term::numresidue;
CTR term::numint;
CTR term::numzeroint;
CTR term::numsimpleint;
CTR term::numdoubleint;
CTR term::nummult;

void term::initialize(int n, int topdeg, bool tracing, 
		      bool factoring, bool zerocheck,
		      METHOD method) {
    U.initialize(topdeg,4);
    Tracing = tracing;
    Factoring = factoring;
    ZeroCheck = zerocheck;
    Method = method;
    TopCoefficient = Method == M_TOPCOEFF;
    S.setup(n,Factoring,TopCoefficient,ZeroCheck);
}