using namespace std;
#include <unistd.h>
#include <getopt.h>

#include "bignum.h"
#include "count.h"
#include "calc.h"
#include "poly.h"
#include "algebra.h"
#include "gen.h"
#include "upoly.h"
#include "version.h"

#include <map>
#include <stack>

int Listing = 0;

void initterm (term& T, int c, int n, int * m) {
    int k;
    term V;

    T = c;

    for(k=1;k<=n;k++) {
	V = zfactor(k,-1,-1);
	T *= V;
    }
    for(k=1;k<=n;k++) {
	for(int j=1;j<=n;j++) 
	{
	    if(j != k) {
		zzfactor zzf = zzfactor(k,j,-m[k]);
		T *= zzf.normalorder();
		V = zzf;
		T *= V;
	    }
	}
	V = zfactor(k,(n-1)*m[k],m[k]);
	T *= V;
    }
    if (Listing) {
	cout << "term (" << m[1];
	for (int p=2;p<=n;p++)
	    cout << " " << m[p];
	cout << "): " << T << "\n";
    }
}

void addterm(int n, int * m, int c, int j, int t, expr& E) {
    term T;
    int k;
    if (j == n) {
	k = n-t;
	m[n] = k;
	initterm(T,c,n,m);
	E += T;
    }
    else {
	int s = t;
	for (k=0; k<= n-t; k++,s++) {
	    m[j] = k;
	    if (k>0)
		c = (c*(n-s+1))/k;
	    addterm(n,m,c,j+1,s,E);
	}
    }
}

// strict = 1 for strict upper triangular, 0 for weak upper triangular
void addupperterm(int n, int * m, int c, int j, int t, int strict, expr& E) {
    term T;
    int k;
    if (j == n) {
	k = n-t;
	m[n] = k;
	initterm(T,c,n,m);
	E += T;
    }
    else {
	int s = t;
	for (k=0; k<= n-t; k++,s++) {
	    m[j] = k;
	    if (k>0)
		c = (c*(n-s+1))/k;
	    if (s>=j+strict)
		addupperterm(n,m,c,j+1,s,strict,E);
	}
    }
}

generator* makeallsum (int n) {
    expr E;
    int * m = new int[n+1];
    addterm(n,m,1,1,0,E);
    delete[] m;
    expgen* GE = new expgen;
    GE->initialize(E);
    return GE;
}

// strict = 1 for strict upper triangular, 0 for weak upper triangular
generator* makeuppersum (int n, int strict) {
    expr E;
    int * m = new int[n+1];
    addupperterm(n,m,1,1,0,strict,E);
    delete[] m;
    expgen* GE = new expgen;
    GE->initialize(E);
    return GE;
}

generator* maketransprod(int c, int m, int n, int* a, int* b) {
    vector<expr> V(n);
    term T;
    zfactor Z;
    zzfactor ZZ;
    expr F;
    for(int k=1; k<=n; k++) {
	F = 0;
	T = 1;
	for(int j = 1; j <= m; j++) {
	    T = zfactor(j,b[k]+m-1,0);
	    for(int i = 1; i <= m; i++) {
		if (i != j) {
		    ZZ = zzfactor(j,i,-1);
		    T *= ZZ.normalorder();
		    T *= ZZ;
		}
	    }
	    F += T;
	}
	V[k-1] = F;
    }

    T = c;
    for(int j=1; j<=m; j++) {
	Z = zfactor(j,-a[j]-1,0);
	T *= Z;
    }

    if (Listing) {
	cout << "Transportation (";
	for (int k = 1; k<= m; k++)
	    cout << a[k] << (k < m ? " " : ") (");
	for (int k = 1; k<= n; k++)
	    cout << b[k] << (k < n ? " " : ")\n");
	cout << T << "\n";
	for(int k=1; k<=n; k++) 
	    cout << "* (\n" << V[k] << ")\n";
    }

    V[0] *= T;
    prodgen* GP = new prodgen;
    GP->initialize(V);
    return GP;
}

int usage(char * prog) {
    cerr << "usage: " << prog << 
	"[options] [path | margins]\n" <<
	"    options: -T -b -n -l -f -z -p -m meth -u n -w n -a n -t -q -h -v\n" <<
       	"             -e range -s num\n" <<
	"    path:    m1 m2 ... mn\n" <<
	"    margins: a1 a2 ... am : b1 b2 ... bn\n";
    return 1;
}

int errmsg(char * mess) {
    cerr << "Error: " << mess << "\n";
    return 2;
}

int main(int argc, char ** argv) {
    int N;
    bool doupper = false;
    bool doall = false;
    bool Polynomial = false;
    int smultiplier = 1;
    int strict = 1;
    bool Transportation = false;
    bool Tracing = false;
    bool Factoring = false;
    bool ZeroCheck = false;
    bool ShowBinomial = false;
    bool ShowNumeric = false;
    bool Quiet = false;
    char* eval_range = 0;
    int eval_start = 1;
    int eval_end = 0;

    CTR numtopc;
    CTR numlowdeg;
    ctr maxstack;
    ctr maxscomp;

    int c;
    char methch = 'b';
    METHOD meth;

    N = 0;
    while ((c = getopt(argc,argv,"u:w:a:fzpbnlm:tqvhe:s:T")) >= 0)
	switch(c)
	{
	case 'u': doupper = true; strict = 1; N = atoi(optarg); break;
	case 'w': doupper = true; strict = 0; N = atoi(optarg); break;
	case 'a': doall = true; N = atoi(optarg); break;
	case 'f': Factoring = true; break;
	case 'z': ZeroCheck = true; break;
	case 'p': Polynomial = true; break;
	case 'b': ShowBinomial = true; break;
	case 'n': ShowNumeric = true; break;
	case 'l': Listing++; break;
	case 'm': methch = *optarg; break;
	case 't': Tracing = true; break;
	case 'q': Quiet = true; break;
	case 'v': cout << VERSION << "\n"; return 0;
	case 'e': eval_range = optarg; break;
	case 's': smultiplier = atoi(optarg); break;
	case 'T': Transportation = true; break;
	case 'h':
	default:
	    return usage(argv[0]);
	}
    
    switch (methch) {
    case 'e': meth = M_EXT; break;
    case 'i': meth = M_INT; break;
    case 'q': meth = M_QUERY; break;
    case 'b': meth = M_BRANCHCOUNT; break;
    case 't': meth = M_TOPCOEFF; break;
    default:
	return usage(argv[0]);
    }

    if (Transportation) {
	Polynomial = false;
    }

    if (N <= 0) {
	if (optind >= argc)
	    return errmsg("nothing to do");
    }

    if (eval_range) {
	eval_start = atoi(eval_range);
	for(;*eval_range && *eval_range!=':'; eval_range++);
	if (*eval_range!=':')
	    return errmsg("evaluation range must be in the form start:stop");
	eval_end = atoi(eval_range+1);
    }

    stack<generator*> GS;
    uint GSheight = 0;

    if (doupper) {
	GS.push(makeuppersum(N,strict));
    }
    else if (doall) {
	GS.push(makeallsum(N));
    }
    else if (Transportation) {
	int n = 0;
	int m = 0;
	int tota = 0;
	int totb = 0;
	bool bstart = false;
	for(int k=optind;k<argc;k++) {
	    if (strcmp(argv[k],":") == 0 || strcmp(argv[k],";") == 0) {
		if (bstart)
		    return errmsg("there must be exactly one ';' character "
				  "separating the margin arrays");
		bstart = true;
	    }
	    else {
		int j = atoi(argv[k]);
		if (j < 0) 
		    return errmsg("margin entry not a non-negative integer");
		if (bstart) {
		    totb += j;
		    n++;
		}
		else {
		    tota += j;
		    m++;
		}
	    }
	}
	if (n == 0 || m == 0)
	    return errmsg("must specify both margin arrays");
	if (tota != totb)
	    return errmsg("margin arrays must have same sum");
	int * a = new int[m+1];
	int * b = new int[n+1];

	for(int k = 1; k <= m; k++)
	    a[k] = atoi(argv[optind + k - 1]);
	for(int k = 1; k<= n; k++)
	    b[k] = atoi(argv[optind + m + k]);
	GS.push(maketransprod(smultiplier,m,n,a,b));
	N = m;

	delete[] a;
	delete[] b;
    }
    else {
	int n = 0;
	for(int k=optind;k<argc;k++) {
	    int j = atoi(argv[k]);
	    if (j < 0) 
		return errmsg("a partition may not have negative terms");
	    n += j;
	}
	if (n != argc - optind)
	    return errmsg("terms do not form a partition");
	int * a = new int[n + 1];
	for(int k=optind;k<argc;k++)
	    a[k+1-optind] = atoi(argv[k]);
	term T;
	initterm(T,smultiplier,n,a);
	delete[] a;
	expr E;
	E = T;
	expgen* GE = new expgen;
	GE->initialize(E);
	GS.push(GE);
	N = n;
    }

    if (Listing > 1)
	return 0;

    maxstack.max(++GSheight);
    
    int TopDeg = (N-1)*(N-1);
    C.precalc(N,TopDeg);

    term::initialize(N,TopDeg,Tracing,Factoring,ZeroCheck,meth);

    poly P(TopDeg);
    
    INT DENOM;
    DENOM = C.Factorial(TopDeg);
    if (!P.fixdenom(DENOM))
	return errmsg("can\'t normalize polynomial P");

    term t;
    tterm tt;

    typedef map<string,tterm> ttMap;
    ttMap M;
    typedef ttMap::const_iterator CTTMI;
	
    RAT topc;
    topc = 0;
    RAT tc;

    INT TRVAL;
    INT TMP;
    TRVAL = 0;

    while (!GS.empty()) {
	if (GS.top()->empty()) {
	    delete GS.top();
	    GS.pop();
	    GSheight--;
	    continue;
	}
	GS.top()->generate(t);
	if (t.constant()) {
	    maxscomp.maxabs(t.scomp());
	    if (Transportation)
		TRVAL += t.evalcoeff(TMP);
	    else if (ShowNumeric) {
		tt = t;
		M[tt.varstr()] += tt;
	    }
	    else {
		numtopc++;
		if (t.topdeg() == TopDeg) {
		    t.topcoeff(tc);
		    if (Tracing) 
			cerr << "topcoeff of " << t << " is " << 
			    tc << "\n";
		    topc += tc;
		}
		else {
		    if (Tracing) 
			cerr << "topcoeff of " << t << 
			    " is of degree " << t.topdeg() << "\n";
		    numlowdeg++;
		}
	    }
	    if (ShowNumeric)
		continue;
	    if (ShowBinomial)
		cout << (t.scomp() > 0 ? "+" : "") << t << "\n";
	    if (Polynomial) {
		poly& T = t.tpoly();
		assert(T.fixdenom(DENOM));
		P += T;
		if (ShowBinomial)
		    cout << "=== " << T << "\n";
		}
	}
	else {
	    generator* G = t.integrate1();
	    if(G) {
		GS.push(G);
		maxstack.max(++GSheight);
	    }
	}
    }  

    if (Transportation)
	cout << "\nEvaluation = " << TRVAL << "\n";
    if (ShowNumeric) {
	topc = 0;
	cout << "\nAfter collecting terms:\n\n";
	for (CTTMI im = M.begin(); im != M.end(); ++im)
	    if (!im->second.zero()) {
		cout << im->second << "\n";
		if (Polynomial) {
		    poly& T = im->second.tpoly();
		    assert(T.fixdenom(DENOM));
		    P += T;
		    cout << "=== " << T << "\n";
		}
		numtopc++;
		if (im->second.topdeg() == TopDeg) {
		    im->second.topcoeff(tc);
		    if (Tracing) 
			cerr << "topcoeff of " << im->second << " is " << 
			    tc << "\n";
		    topc += tc;
		}
		else {
		    if (Tracing) 
			cerr << "topcoeff of " << im->second << 
			    " is of degree " << im->second.topdeg() << "\n";
		    numlowdeg++;
		}
	    }
    }
    if(Polynomial)
	cout << "\nFull polynomial = \n" << P << "\n";
    if(!Transportation)
	cout << "\nTop coefficient = " << topc << "\n\n";

    if(Polynomial && eval_start <= eval_end) {
	cout << "Polynomial values:\n";
	RAT v;
	for(int k=eval_start; k<=eval_end; k++)
	    cout << "    P(" << k << ") = " << P.eval(v,k) << "\n";
    }

    if(Quiet)
	return 0;
    cout << "\nVersion " << VERSION << "\n";
    cout << "\nnumber of integrals evaluated:  " << term::numint << "\n";
    cout <<   "    number of zero integrals:   " << term::numzeroint << "\n";
    cout <<   "    number of simple integrals: " << term::numsimpleint << "\n";
    cout <<   "    number of double integrals: " << term::numdoubleint << "\n";
    cout <<   "    number of factorizations:   " << term::nummult << "\n";
    cout <<   "number of residues calculated:  " << term::numresidue << "\n";
    cout <<   "number of top coeff terms:      " << numtopc << "\n";
    cout <<   "    number of low degree:       " << numlowdeg << "\n";
    if (Polynomial) {
	cout << "\n";
	cout << "polynomial cache size = " << upoly::numsetcache << "\n";
	cout << "           cache hits = " << upoly::numgetcache << "\n";
	cout << "\n";
    }

//     cout << "\nmax factorial argument = " << C.maxfact_n << "\n";
//     cout << "max binomial argument = " << C.maxibin_n << "\n";
//     cout << "max binomial value = " << C.maxibin_val << "\n";
    cout << "max nmc arguments = " << C.maxnmc_r << ", " << C.maxnmc_m << "\n";
    cout << "max nmc value = " << C.maxnmcval << "\n";
    cout << "max scalar factor = " << maxscomp << "\n";
    cout << "max GS stack height = " << maxstack << "\n";

    return 0;
}