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

proc remove (list l, int x) {

list rl;

for (int i=1; i<=size(l); i++) {

	if (l[i][1]<>x) {
		rl = insert(rl,l[i]);
	}
}

rl = reverse (rl);

return(rl);

}

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

proc listify (ideal I) {

list l;

if (I==0) {
	return (list(0));
};

for (int i=1; i<=size(I); i++) {
	l[i]=I[i];
}
return(l);
};

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

proc flatlist (list l) {

list resl;

for (int i=1; i <= size(l); i++ ) {
	
	resl = resl + l[i];

}

return (resl);

};

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

proc allVarProd () {

poly p = 1;
for (int i=1; i <= nvars(basering); i++ ) {
	p= p * var(i);
};

return(p);
};
			

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

proc fplace (poly f, list l) {

//print ("fplace: = ");
//print(l);

for (int i=1; i <= size(l); i++ ) {

	if (l[i]==f) {
		return(i);
	};

}

print("fplace: not found");
return(-99);

}

			
// matrix Pjkmat = polyMultMat(si1sum[j],sisum[k],mri[j,k]);


proc polyMultMat(list monc, list monr, poly g) {

int ncol = size(monc);
int nrow = size(monr);

matrix Pjk[nrow][ncol];

//print("polyMultMat:monr = ");
//print(monr);

poly mon, h;
matrix mh;
int p;
int s;

for (int i=1; i<=ncol; i++) {

	mon = monc[i];
	h = mon * g;

	if (not (h==0)) {
		mh = coef (h, allVarProd());

	//	print("mh=");
	//	print(mh);

		for (s=1; s <= ncols(mh); s++ ) {

			p = fplace (mh[1,s], monr);
			Pjk[p,i] = mh[2,s];


		};
	}

};

return(Pjk);

};

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

proc mmaxideal (int k) {


list lm;

if (k < 0) {
	lm = list(0);
} else  {
	lm = listify(maxideal(k));
};

return(lm);
};


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

proc repeatn (int val, int n) {

list l;
for (int i=1; i <= n; i++) {
	l[i]=val;
}
return(l);
}

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

proc myfac (int n) {

number i = n;

number p = 1;

while (i > 1) {
	p = p * i;
	i = i - 1;
};

return(p); 

};



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

proc dimOXd (int n, int d) {

//print("n/d=");print(n);print(d);

if (d < 0) {
	return ( 0 );
};

int zz = n - 1 + d;
int dd = n - 1;


//print ("zz/dd=");
//print (zz); print (dd);
string rx=binomial(zz,dd);

execute("int rr= " + rx );

kill rx;

return (rr);

}

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

proc degDiff (list lldeg, matrix mm, int ncol ) {

int d = 0;
/*
print("mm=");
print(mm);
print (lldeg);
*/
for (int i=1; i <= nrows(mm); i++ ) {

	d = lldeg[i];
	
	if (mm[i,ncol] <> 0) {
		return ( d + deg(mm[i,ncol]) );
	};

};

return(d);

};

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

proc sheafCohJB ( module M, int dlow, int dhigh ) {

int n=nvars(basering);

int dd = dhigh - dlow + 1;

intmat cohMat[n][dd];

list resl;
int kk;

for (int d=dlow; d <= dhigh; d++ ) {

	resl = sheafCohJBx (M, d);
	for ( kk= 1; kk <= size(resl); kk++) {
		cohMat[n-kk+1,d-dlow+1] = resl[kk];
	};
};

return (cohMat);

};


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

proc sheafCohJBx (module M, int d)
{

def b=basering;

int n = nvars(basering);
int dimP = n - 1;


resolution ri = mres(M,0);

list rilis;

list hdvl;

list hdvli;

int imax;

intmat bri = betti(ri);

int jj;

for (int i=1; i <= ncols(bri); i++ ) {

	//print ("i=");
	//print(i);

	hdvli = list();
	for (jj=1; jj <= nrows(bri); jj++) {
		if (bri[jj,i]<>0) {
			hdvli=hdvli+repeatn(i+jj-2,bri[jj,i]);
		};
	};
		
	hdvl[i] = hdvli;

	if (i-1 >= 1) {
		rilis[i-1] = matrix(ri[i-1]);
	};
};

imax = i - 1;
int ii;

for ( i=i; i <= n+1; i++ ) {

	hdvli = list();
	
	// print ("nrows/ncols=");print(nrows(matrix(ri[i-1])));print(ncols(matrix(ri[i-1])));

	rilis[i-1] = matrix(ri[i-1]);
	
	//kill mz00;

	for (ii=1; ii <= ncols(rilis[i-1]); ii++ ) {
	
		hdvli = insert (hdvli, degDiff(hdvl[i-1],rilis[i-1],ii));
	
	}

	hdvl[i] = reverse ( hdvli );

	//kill rlisa;

};

for ( i=i; i <= n+3; i++ ) {

	hdvli = list();
	
	// print ("nrows/ncols=");print(nrows(matrix(ri[i-1])));print(ncols(matrix(ri[i-1])));


	matrix mz00[nrows(rilis[i-2])][1]=0;

	matrix rlisa[ncols(rilis[i-2])][1]=matrix(modulo(rilis[i-2],mz00));

	rilis[i-1] = rlisa;
	
	kill mz00;

	for (ii=1; ii <= ncols(rilis[i-1]); ii++ ) {
	
		hdvli = insert (hdvli, degDiff(hdvl[i-1],rilis[i-1],ii));
	
	}

	hdvl[i] = reverse ( hdvli );

	kill rlisa;

};

int nrilis=size(rilis);

/*
print ("hdvl=");

print (hdvl);
print( rilis );
*/

list matrlis;
list bigAlis;

int j, k;
int rki1, rki;
int nrowsPjkmat;
int ncolsPjkmat;
int bigarow;
int bigacol;

list si1sum;
list sisum;

list i1mons;
list imons;


int nni1;
int nni;

matrix Pjkmat;
// matrix mri;

list exZeroLis;

for (i=1; i <= nrilis+1; i++) {
	exZeroLis[i]=0;
};


for (i=1; i<= nrilis; i++ ) {

	matrix mri = rilis[i];
	
	//print("i=");print(i);
	
	//print ("mri= ");print(mri );
	
	
	rki1 = nrows(mri);
	rki = ncols(mri);
/*
	print ("row/col=");
	print (rki1);
	print(rki);
	
	print ("hdvli=");
	print(hdvl[i]);
	print(hdvl[i+1]);
*/
	si1sum = list();
	sisum = list();
	
	i1mons = list();
	imons = list();
	
	//print ("j's = "); 
	
	for (j=1; j <= rki1; j++ ) {
		
		if (hdvl[i][j]-n-d < 0) {
			exZeroLis[i]=exZeroLis[i]+1;
		};
		
		si1sum[j]=mmaxideal(hdvl[i][j]-n-d);

		//print(si1sum[j]);

		i1mons = insert(i1mons,si1sum[j]);

	};
	
	//print ("k's = ");
	
	for (k=1; k <= rki; k++ ) {

		if ((i == (nrilis-1)) and (hdvl[i+1][k]-n-d < 0)) {
			// print ("now/i = ");print(i);
			exZeroLis[i+1]=exZeroLis[i+1]+1;
		};
	
		sisum[k]=mmaxideal(hdvl[i+1][k]-n-d);
		
		//print(sisum[k]);
		
		imons = insert(imons,sisum[k]);
	
	};
	
	i1mons = reverse (i1mons);
	imons = reverse (imons);
/*
	print ("i1/imons = ");
	print (i1mons);
	print (imons);
*/	

	nni1 = size (flatlist(i1mons));
	nni = size(flatlist(imons));
			
	matrix biga[nni][nni1];
	
	bigacol = 1;
	
	//print ("nni1/nni=" + string(nni1)+ " " + string(nni));
	
	//print ("mri= ");
	//matprint ( mri );
		
	for (j=1; j <= rki1; j++ ) {
	
		bigarow = 1;
	
		for (k=1; k <= rki; k++ ) {
	
			Pjkmat = polyMultMat(si1sum[j],sisum[k],mri[j,k]);

			//print ("Pjkmat = " );
			//print ( Pjkmat );

			nrowsPjkmat = nrows(Pjkmat);
			ncolsPjkmat = ncols(Pjkmat);

			
			biga[bigarow..bigarow+nrowsPjkmat-1,bigacol..bigacol+ncolsPjkmat-1] = Pjkmat;
			
			bigarow = bigarow + size(sisum[k]);

		};
		
		bigacol = bigacol + size(si1sum[j]);
	
	};


	bigAlis = insert(bigAlis, biga);
	
	kill biga;
	kill mri;

};

//print(exZeroLis);

bigAlis = reverse(bigAlis);


for (i=1; i <= size(bigAlis); i++ ) {
	bigAlis[i]=transpose(bigAlis[i]);
};

int nmax = size(bigAlis);
/*
matrix mz00[ncols(bigAlis[nmax])][1]=0;

bigAlis[nmax+1]=mz00;
*/

int rrk = mat_rk(bigAlis[3]);
/*
print("nrows=");print(nrows(bigAlis[3]));
print("ncols=");print(ncols(bigAlis[3]));
print("rrk=");print(rrk);
*/

//print(size(bigAlis));

list reslis;

matrix mz0[1][nrows(bigAlis[1])]=0;

reslis=insert (reslis, matrix(homology(bigAlis[1],mz0,0,0)));

matrix homol;

//print ("*** Computing Homology");

// print ("size(bigAlis) = ");print(size(bigAlis));

for (i=1; i <= size(bigAlis)-1; i++) {
	
	//print ("i=");print(i);
	homol = matrix(homology(bigAlis[i+1],bigAlis[i],0,0));
	//print ("homol=");print(homol);
	reslis = insert(reslis, homol);
};

reslis = reverse ( reslis );

matrix rli;

//print ("size(reslis)=");print(size(reslis));

int eZL = 0;

for (i=1; i <= size(reslis); i++ ) {
	//print ("dim: i="); print(i);
	
	rli = reslis[i];
	
	eZL = exZeroLis[i];
	
	reslis[i] = nrows(rli)-mat_rk(rli)-eZL;
};

//print ("*** Computing Zero Cohomology");

int hnL = reslis[n];

rli = bigAlis[size(bigAlis)];

int hn1L = reslis[n+1];

//print ("hnL = ");print(hnL);
//print ("hn1L = "); print(hn1L);
//print ("reslis = "); print(reslis);

int s = 0;
int s1 = 0;
int eps = 1;
int kk;
int slast;

for (int ll=1; ll <= imax; ll++) {

	s = 0;
	for (kk=1; kk <= size(hdvl[ll]); kk++ ) {
		s = s + dimOXd ( n, -hdvl[ll][kk] + d );
	};
	slast = s;
	s1 = s1+ eps * s;
	eps = -eps;
};

/*
if (matrix(ri[n])==0) {
	s1 = s1 + eps * slast;
}
*/

//print ("s1=");print(s1);

int hom0dim = s1 + hnL - hn1L;

reslis = reslis[1..n-1];

reslis = reslis + list(hom0dim);

reslis = reverse ( reslis );


//print("reslis = " + string(reslis));

return ( reslis );

};



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

LIB "sheafcoh.lib";


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

proc testJB () {

string str="ring r=0,(x,y,z,w),(dp(4))";

execute(str);

export(r);

ideal i0=0;
ideal i1=x3y+w4+z2x2;
ideal i2=x2+y2,xy+zw;

module m0 = i0;
module m1 = i1;
module m2 = i2;

homog(m0);
homog(m1);
homog(m2);

module merr=module(
15826x2y*gen(1)+12074y2z*gen(1)+10614xz2*gen(1)+19435y2w*gen(1)+23389xw2*gen(1)+
21369yw2*gen(1)+6251w3*gen(1),
22843y2z*gen(1)+10165xz2*gen(1)+12008x2w*gen(1)+29369xyw*gen(1)+18277yw2*gen(1)+13880w3*gen(1),
3483xzw*gen(1)+14644w3*gen(1));

export(i0,i1,i2,m0,m1,m2,merr);

keepring r;

};

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

proc compChi (int n, int d) {

	int erg = dimOXd ( n, d ) + (-1)^(n-1)* dimOXd ( n, -d - n );

	return (erg);
	
};


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

proc correctZeroLine (intmat cohtab, module mtst, int dlow ) {

int n = nvars(basering);

resolution rtst = mres(mtst,0);

intmat btst = betti (rtst);

int zrow = nrows(cohtab);

int bi, bj;

int eps = 1;
int s, s1;

int dakt;

for (int acol=1; acol <= ncols(cohtab); acol++ ) {

	dakt = dlow + acol - 1;

	s = 0;
	eps = 1;

	for (bj = 1; bj <= ncols (btst); bj++ ) {

		s1 = 0;

		for (bi =1; bi <= nrows (btst); bi++ ) {

			s1 = s1 + btst[bi,bj] * compChi (n, -(bi+bj-2)+dakt );


		};

		s = s + eps * s1;

		eps = -eps;

	};

	eps = 1;

	for (bi=nrows(cohtab)-1; bi >= 1; bi-- ) {
	
		s = s + cohtab[bi,acol]*eps;
		
		eps=-eps;
	};
	cohtab[nrows(cohtab),acol]=s;

};




return (cohtab);

};


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


proc verifyJB ( int n, int anzTests ) {

string vs = "(x(1.."+string(n)+"))";

string str="ring r=0,"+vs+",dp";


execute(str);

setring r;

keepring(r);

intmat cohtabjb;
intmat cohtab;

ideal idtst;
module modtst;

int errs = 0;
list errlis;

int ngenId;
int degId;

for (int tcnt=1; tcnt <= anzTests; tcnt++ ) {

	ngenId=random(1,n-1);
	degId = random(1, 4);

	idtst = sparseid(ngenId, degId,degId, 75, 5);

	if (tcnt==1) {
		idtst=0;
	};

	modtst=module(idtst);
	if (homog(modtst) <> 1) {
		print ("homogeneity failed");
	};
	
	cohtab = sheafCoh(modtst,-5,5);
	
	cohtab=correctZeroLine(cohtab,modtst, -5);
	
	print("cohtab="); print(cohtab);
	
	cohtabjb = sheafCohJB(modtst,-5,5);
	
	
	
	if (not(cohtab==cohtabjb)) {
		errs++;
		errlis=insert(errlis,modtst);
		print ("Test: ** Error.");
		return (modtst);	
	} else {
		print ("Test: Ok.");
	};
		
};

return (errs, errlis);
};

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

proc matprint (matrix m) {

int nrow = nrows(m);
int ncol = ncols(m);

int j;

for (int i=1; i<=nrow; i++ ) {
	for (j=1; j<=ncol; j++ ) {
		print("i/j="+string(i)+" "+string(j));
		print(m[i,j]);	
	};
};

}