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

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

evaluate maps the variables z_k that are present in the term t
to internal variables w_j, and sets up various arrays:

IK[j] is a Z iterator for the variable w_j
KV[k] is the w index for the variable z_k
Q[k1][k2] is the exponent of the (w_k1 - w_k2) factor
NZZ[k] is the number of ZZ factors involving w_k
DZZ[k] is the degree of the ZZ factors involving w_k
E contains the evaluation data for each group

a group is a consecutive list of variables that can be grouped for 
integration.  groups will have size 1 if the Factoring flag is false

evaluate depends on the TopCoefficient flag only through the maxzdiff
call, which computes the maximum number of times a zfactor should be
differentiated during a residue calculation (see description of 
zfactor::maxzdiff()

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

void statterm::evaluate(const term &t) {
    nvar = 0;
    for(int v = 1; v<=N; v++)
	KV[v] = 0;
    for(CZI iz = t.Z.begin(); iz != t.Z.end(); ++iz) {
	KV[iz->var()] = ++nvar;
	IK[nvar] = iz;
    }
    for(int k1 = 1; k1<=nvar; k1++) {
	NZZ[k1] = 0;
	DZZ[k1] = 0;
	for(int k2 = 1; k2<=nvar; k2++)
	    Q[k1][k2] = 0;
    }
    for(CZZI izz = t.ZZ.begin(); izz != t.ZZ.end(); ++izz) {
	int k1 = KV[izz->var(1)];
	int k2 = KV[izz->var(2)];
	int e = izz->exp();
	Q[k1][k2] = e;
	Q[k2][k1] = e;
	NZZ[k1]++;
	NZZ[k2]++;
	DZZ[k1] += e;
	DZZ[k2] += e;
    }

/********************
  tries to prune zero integrals as early as possible; this can be done when
  total degree in the first r variables is less than -r, since further 
  external integrations among the first r variables will raise this by 1 --
  so eventually all integrals arising from this will have degree < -2 in the 
  first integral

  this is only a substantial improvement on large integrals (n=10, I hope)
********************/
  
    if (ZeroCheck) {
	lin T;
	zeroz = 0;
	for(int r=1; r<nvar; r++) {
	    T += IK[r]->exp();
	    int qsum = 0;
	    for(int k=r+1; k<=nvar; k++)
		qsum += Q[r][k];
	    T += qsum;
	    if (T < -r) {
		zeroz = r;
		return;
	    }
	}   
    }

    nblock = 0;
    lin e;
    for(int k = 1; k<= nvar; ) {
	e = IK[k]->exp();
	int p;

	p = IK[k]->maxzdiff(INT_MAX,TopCoefficient);
	E[++nblock].start = k;

// int_ok is just term::zero_res_at_zero()
	E[nblock].int_ok = !(e < 0);
	if(E[nblock].int_ok) {
	    int tot = 0;
	    for(int i = nvar; i>k; i--)
		if(Q[i][k])
		    tot += C.nmc(-Q[i][k]-1,NZZ[k],p);
	    E[nblock].int_eval = tot;
	}
	else
	    E[nblock].int_eval = -1;

// ext_ok is just term::zero_res_at_infinity()
	E[nblock].ext_ok = e.tpart() < 0 || e.tpart() == 0 && 
	    e.ipart() + DZZ[k] < -1;
	if(E[nblock].ext_ok) {
	    int tot = 0;
	    for(int i = 1; i<k; i++)
		if(Q[k][i])
		    tot += C.nmc(-Q[k][i]-1,NZZ[k],p);
	    E[nblock].ext_eval = tot;
	}
	else
	    E[nblock].ext_eval = -1;
	
	int j = k+1;
	if (Factoring) {
	    for(; j<= nvar; j++) {
		if(!(IK[j]->exp()==e))
		    break;
		bool OK = true;
		for(int i = 1; i <= nvar; i++)
		    if(Q[i][j] != Q[i][k]) {
			OK = false;
			break;
		    }
		if (!OK)
		    break;
	    }
	}

	E[nblock].mult = j-k;
	k=j;
    }
    besteval = INT_MAX;
    bestdir = ' ';
    int beststart = 0;
    bestmult = 0;
    for(int k=1; k<=nblock; k++) {
	if(E[k].mult > bestmult && (E[k].ext_ok||E[k].int_ok)) {
	    bestmult = E[k].mult;
	    beststart = k;
	    if (E[k].ext_ok) {
		bestdir = 'e';
		besteval = E[k].ext_eval;
	    }
	    else {
		bestdir = 'i';
		besteval = E[k].int_eval;
	    }
	    continue;
	}
	if (E[k].mult < bestmult)
	    continue;
	if (E[k].ext_ok && E[k].ext_eval < besteval) {
	    beststart = k;
	    bestdir = 'e';
	    besteval = E[k].ext_eval;
	}
	if (E[k].int_ok && E[k].int_eval < besteval) {
	    beststart = k;
	    bestdir = 'i';
	    besteval = E[k].int_eval;
	}
    }
    bestz = IK[E[beststart].start];
}

bool statterm::ok_int(int v, char d, int m, CZI& qz) {
    if(v<1 || v>N)
	return false;
    int k = KV[v];
    if(!k)
	return false;
    qz = IK[k];
    if(d != 'e' && d != 'i')
	return false;
    if(m<1)
	return false;
    int n;
    for(n=nblock; k<E[n].start; n--);
    if(d == 'i' && !E[n].int_ok)
	return false;
    if(d == 'e' && !E[n].ext_ok)
	return false;
    if(k+m-1 > E[n].start+E[n].mult-1)
	return false;
    return true;
}

void statterm::setup(int n, 
		     bool factoring, bool topcoefficient, bool zerocheck) {
    Factoring = factoring;
    TopCoefficient = topcoefficient;
    ZeroCheck = zerocheck;
    zeroz = 0;

    N = n;
    IK = new CZI [N+1];
    KV = new int [N+1];
    Q = new int* [N+1];
    for(int k1 = 1; k1 <=N; k1++) {
	Q[k1] = new int [N+1];
	Q[k1][k1] = 0;
    }
    NZZ = new int [N+1];
    DZZ = new int [N+1];
    E = new evaldata[N+1];
}

ostream& statterm::dumpeval(ostream& f) const {
    if(ZeroCheck && zeroz) {
	f << "zero integral based on first " << zeroz <<" vars\n";
	return f;
    }
    for(int k=1; k<=nblock; k++)
	f << "  start z_" << IK[E[k].start]->var() <<
	    "; end z_" << IK[E[k].start]->var()+E[k].mult-1 <<
	    " (" << E[k].mult << ")" <<
	    (E[k].int_ok ? " :i" : " : ") <<
	    (E[k].ext_ok ? "e: " : " : ") <<
	    "int eval = " << E[k].int_eval <<
	    "; ext eval = " << E[k].ext_eval << "\n";
    f << "  best: " << bestz->var() << bestdir <<
	" (" << bestmult << "), eval = " << besteval << "\n";
    return f;
}    

ostream& operator<<(ostream& f, const statterm& S) {
    f << "nvar = " << S.nvar << "\n";
    for(int v=1; v<=S.N; v++)
	if(S.KV[v])
	    f << "  z_" << v << " -> w_" << S.KV[v] << "\n";
    for(int k=1; k<=S.nvar; k++)
	f << "  w_" << k << "^" << S.IK[k]->exp() << "; ";
    f << "\n";
    for(int k1 = 1; k1<=S.nvar; k1++)
	for(int k2 = 1; k2<k1; k2++)
	    if(S.Q[k1][k2])
		f << "(w_" << k1 << " - w_" << k2 << ")^" <<
		    S.Q[k1][k2] << "; ";
    f << "\n";
    for(int k=1; k<=S.nvar; k++)
	f << "  w_" << k << " zzfactors = " << S.NZZ[k] << 
	    "; zzdegree = " << S.DZZ[k] << "\n";
    f << "nblock = " << S.nblock << "\n";
    return S.dumpeval(f);
}