newPackage(
	"MultiProjChar",
	Version => "1.2", 
    	Date => "April 24, 2015",
    	Authors => {{Name => "Martin Helmer", 
		  Email => "martin.helmer2@gmail.com", 
		  HomePage => "http://publish.uwo.ca/~mhelmer2/"}},
    	Headline => "Computes CSM classes, Segre classes and the Euler Char. in Multi-Projective Space",
    	DebuggingMode => false,
	Reload => true
    	)

needsPackage "NormalToricVarieties";

--Exported functions/variables
export{
   Segre,
   CSM,
   Euler,
   ChowRing,
   Output,
   ListForm,
   ChowRingElement,
   Method,
   InclusionExclusion,
   DirectCompleteInt

}


ChowRing=method(TypicalValue=>QuotientRing);
ChowRing NormalToricVariety:=TorVar->(
    assert isSimplicial TorVar;
    --First build Chow ring, need Stanley-Reisner Ideal (SR) and the ideal
    -- generated by the linear relations of the rays (J)
    --See Cox, Little, Schenck Th. 12.5.3 and comments after proof
    R:=ring(TorVar);
    --For simplical toric var. Lemma 3.5 of Euler characteristic of coherent sheaves on simplicial torics via the Stanley-Reisner ring
    -- (and probally other sources) tell us that the SR ideal is the Alexander
    --dual of the toric irrelevant ideal 
    SR:=dual monomialIdeal TorVar;
    F:=fan TorVar;
    Fd:=dim(F);
    --Build ideal generated by linear relations of the rays
    Jl:={};
    for j from 0 to dim(F)-1 do(
	Jl=append (Jl,sum(0..((length rays(TorVar))-1), i->(((rays TorVar)_i)_j)*(gens R)_i ));
    	);
    J:=ideal(Jl);
    --Chow ring
    A:=R/(SR+J);
    --Generators (as a ring) of the quotient ring representation of the chow ring corespond to 
    --the divisors associated to the rays in the fan Theorem 12.5.3. Cox, Little, Schenck and 
    --comments above
    return A;
    
    
    );
ChowRing Ring:=R->(
    Rgens:=gens R;
    Rdegs:=degrees R;
    degd:=0;
    eqs:=0;
    ChDegs:=unique Rdegs;
    m:=length ChDegs;
    C:=ZZ[gens R, Degrees=>degrees R, Heft=>heft R];
    K:={};
    inds:={};
    rg:=0;
    ns:={};
    temp:=0;
    for d in ChDegs do(
	temp=0;
	for a in Rdegs do(
	    if d==a then temp=temp+1;
	    );
	ns=append(ns,temp);
	);
    
    for i from 0 to length(ChDegs)-1 do(
	eqs=m->m==ChDegs_i;
	inds=positions(Rdegs,eqs);
	degd=Rgens_inds;
	K=append(K,degd_0^(ns_i));
	for j from 1 to length(degd)-1 do(
	    K=append(K,degd_0-degd_j);
	    );
	
	);
    K=substitute(ideal K,C);
    A:=C/K;
    return A;
);
Euler = method(TypicalValue => RingElement,Options => {Method=>InclusionExclusion});

Euler Ideal:=opts->I->(
    A:=ChowRing(ring I);
    B:=last flatten entries sort basis A;
    csm:=CSM(A,I,Method=>opts.Method);
    EC:=csm_(B);
    <<"The Euler Chacteristic = "<<EC<<endl;
    return EC;
    );
Euler RingElement:=opts->csm->(
    <<"Finding the Euler Chacteristic using input csm class= "<<csm<<endl;
    A:=ring csm;
    B:=last flatten entries sort basis A;
    EC:=csm_(B);
    <<"The Euler Chacteristic = "<<EC<<endl;
    return EC;
    );

CSM = method(TypicalValue => RingElement,Options => {Method=>InclusionExclusion});

CSM (NormalToricVariety, Ideal):= opts->(TorVar,I)->(
    return CSM(ChowRing(TorVar),I,Method=>opts.Method);
    );
CSM Ideal:=opts->I->(
    return CSM(ChowRing(ring I),I,Method=>opts.Method);
    );

CSM (QuotientRing,Ideal):= opts -> (ChRing,I)->(
    R:=ring I;
    A:=ChRing;
    gensI:=flatten entries gens I;
    B:=flatten entries sort basis A;
    ns:=degree last B;
    n:=sum(ns);
    m:=length ns;
    SegY:=0;
    SegV:=0;
    V:=0;
    K:=0;
    J:=0;
    vf:=0;
    csm:=0;
    csm2:=0;
    ChernTR:=product(m,i->(1+(basis(OneAti(m,i),A))_0_0)^(ns_i+1));
    if opts.Method==DirectCompleteInt then(
	<<"Trying direct"<<endl;
	use R;
	codimI:=codim I;
	irelId:=irrell(R);
	J2:=0;
	Z:=0;
	Sing:=0;
	cont:=false;
	if codimI!=length(gensI) then(
	    <<"Input is not a complete intersection, using inclusion/exclusion instead"<<endl;
	    return CSM(ChRing,I);
	    )
	else(
	    --check if input smooth?
	    ssets:=subsets(gensI,length(gensI)-1);
	    for s in ssets do (
	    	J2=saturate(minors(length(s),jacobian ideal(s))+ideal(s),irelId);
	    	if codim(J2)>n then (
	    	    Z=ideal(s);
	    	    Sing=ideal toList(set(gensI)-set(s));
	    	    cont=true;
	    	    break;
	    	    );
	    	);
	    if cont then(
		<<"Assumtions okay"<<endl;
		if length(gensI)>1 then J=saturate(minors(length(gensI),jacobian I)+I,irelId);
		Vlist:={};
		dv:=0;
		for f in gensI do(
		    dv=degree f;
	    	    Vlist=append(Vlist,sum(length dv, i-> dv_i*(basis(OneAti(m,i),A))_0_0));
		    );
		r:=length(gensI);
		CE:=product(r,j->(1+Vlist_j));
		dv=degree Sing_0;
		Vr:=sum(length dv, i-> dv_i*(basis(OneAti(m,i),A))_0_0);
		V2:=product(r,j->Vlist_j);
		if r==1 then return (ChernTR*V2)//CE;
		SegY=Segre(A,J,Output=>ListForm);
		CEotimesL:=sum(0..r,i->(-1)^i*(1+Vr)^(r-i)*sum(select(terms CE, q->sum(degree(q))==i)));
		segstuff:=sum(0..n,i->((-1)^i*SegY_(i))//( (1+Vr)^i) );
		cfj:=(ChernTR//CE)*V2;
		milnor:=(ChernTR//CE)*(-1)^(r+1)*CEotimesL*(segstuff);
		csm=cfj+milnor;
		return csm;
		)
	    else(
		 <<"Input does not satisfy assumtions, using inclusion/exclusion instead"<<endl;
	    	 return CSM(ChRing,I);
		);
	);
    )
    else (
	Isubsets:=delete({},subsets(gensI));
	for f in Isubsets do(
	    K=ideal product(f);
	    J=ideal(delete(0_R,flatten entries jacobian K));
	    vf=degree K_0;
	    V=sum(length vf, i-> vf_i*(basis(OneAti(m,i),A))_0_0);
	    J=J+K;
	    SegY=Segre(A,J,Output=>ListForm);
	    SegV=V//(1+V);
	    csm2=(-1)^(length(f)+1)*(ChernTR*(SegV+sum(n+1,i->sum(n+1-i,j->binomial(n-i,j)*(-V)^j*(-1)^(n-j-i)*SegY_(n-i-j)))) );
	    csm=csm+csm2;
	    );
	<<"The CSM class= "<<csm<<endl;
	return csm;
	);
    
    )



Segre = method(TypicalValue => RingElement,Options => {Output=>ChowRingElement});

Segre (NormalToricVariety, Ideal):= opts->(TorVar,I)->(
    return Segre(ChowRing(TorVar),I,Output=>opts.Output);
    );
Segre Ideal:=opts->I->(
    return Segre(ChowRing(ring I),I,Output=>opts.Output);
    );
Segre (QuotientRing,Ideal):= opts->(ChRing,I)->(
    R:=ring I;
    kk:=coefficientRing R;
    A:=ChRing;
    B:=flatten entries sort basis A;
    ns:=degree last B;
    n:=sum(ns);
    deg1B:={};
    for w in B do if sum(degree(w))==1 then deg1B=append(deg1B,w);
    ChDegs:=unique degrees R;
    m:=length ChDegs;
    seg:=0;
    t1:=symbol t1;
    S:=kk[gens R, t1];
    codimI:=codim(I);
    dimI:=n-codimI;
    gensI:= flatten sort entries gens I;
    degI:= degrees I;
    transDegI:= transpose degI;
    len:= length transDegI;
    maxDegs:= for i from 0 to len-1 list max transDegI_i;
    cOMaxDegs:=sum(length deg1B,i->(basis(OneAti(m,i),A))_0_0*maxDegs_i);
    J:= for i from 1 to n list sum(gensI,g -> g*random(maxDegs-degree(g),R));
    RdList:={};
    GList:={};
    Jd:=0;
    JT:=0;
    c:={};
    v:=0;
    K:=0;
    Ls:=0;
    LA:=0;
    gbWt2:=0;
    tall2:=0;
    Yiota:=0;
    if codimI<=n then(
	GList=for iota from 0 to (codimI-1) list cOMaxDegs^(iota);
	--time(
	for iota from codimI to n do(
	    Jd=substitute(ideal take(J,iota),S);
	    JT=ideal (1-t1*substitute(sum(gensI,g -> g*random(kk)),S));
    	    c={};
    	    for w in B do if sum(degree(w))==iota then c=append(c,w);
    	    Yiota=0;
    	    for w in c do(
	    	Ls=0;
	    	K=0;
	    	v=ns-degree(w);
	    	for i from 0 to length(v)-1 do(
	    	    if v_i!=0 then (
		    	Ls=Ls+sum(v_i,j->ideal(random(OneAti(m,i),R)));
		    	);
	    	    );
		--time(
	    	LA=substitute(sum(length deg1B,j->ideal(1-random(degree(deg1B_j),R))),S);
	    	K=Jd+JT+substitute(Ls,S)+LA;
	    	gbWt2 = groebnerBasis(K, Strategy=>"F4");
            	tall2 = numColumns basis(cokernel leadTerm gbWt2);
		--);
	    	Yiota=Yiota+tall2*w;
		--<<tall2<<", ";
	    	);
    	    GList=append(GList,Yiota);
    	    );
	--);
	--the following preforms the aluffi tensor notation comp    
	--GxOMD:=sum(0..n,i->GList_i//((1+cOMaxDegs)^i));
	temp3:=1;
	GxOMD:=0;
	--time(
	for i from 0 to n do(
	    GxOMD=GxOMD+GList_i*(temp3);
	    temp3=temp3//(1+cOMaxDegs);
	    );
	seg=1-GxOMD//(1+cOMaxDegs);
	--);
        --<<"Done Quotient Ring calcs "<<endl;
	tseg:=terms(seg);
	tot:=0;
    	if opts.Output==ListForm then( 
	    segList:={};
	    for i from 0 to n do(
	    tot=0_A;
	    for f in tseg do(
		if sum(degree(f))==i then(
		    tot=tot+f
		    );
	    );
	    segList=append(segList,tot);
	    );	
	    return segList
	    )
	else (
	    <<"Segre= "<<seg<<endl;
	    return seg
	    );
	)
    else(
	segList=for i from 0 to n list 0_A;
	seg=0_A;
	if opts.Output==ListForm then(
	    return segList;
	    )
	 else (
	     <<"Segre= "<<seg<<endl;
	     return seg;
	     );
	);
);

---------------------------
--Internal functions 
---------------------------

OneAti=(dl,i)->(
    vec:={};
    for j from 0 to dl-1 do(
	if j==i then vec=append(vec,1) else vec=append(vec,0);
	);
    return vec;
    )

irrell=R->(    
    Rgens:=gens R;
    Rdegs:=degrees R;
    bloks:=unique Rdegs;
    irId:=ideal 1_R;
    elList:={};
    for a in bloks do(
	elList={};
	for r in Rgens do(
	    if degree(r)==a then(
		elList=append(elList,r);
		);
	    );
	irId=irId*ideal(elList)
	
	);
    return irId;
    )
----------------------------
--Miscilaneous Examples
--
---------------------------

TEST ///
{*
    restart
    needsPackage "MultiProjChar"
*}
   kk=ZZ/32749;
   R = kk[x_0..x_3,y_0..y_2,z_0,z_1,Degrees=>{{1,0,0},{1,0,0},{1,0,0},{1,0,0},{0,1,0},{0,1,0},{0,1,0},{0,0,1},{0,0,1}},Heft=>{1,1,1}]
   I = ideal(x_0*z_0^2,y_0*z_0+z_0*y_1)
   A=ChowRing(R);
   se=time Segre(A,I)
   time CSM(I)
///

TEST ///
{*
    restart
    needsPackage "MultiProjChar"
*}
    kk=ZZ/32749;
    R=kk[x_0..x_6];
    I=ideal ((x_0*x_3-x_1*x_2)*(x_0-x_6),x_0^2*x_5*x_4-x_4*x_0^2*x_3);
    time Segre(I)
    time CSM(I)
///
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
*}
   kk=ZZ/32749;
   R = kk[x_0..x_2,y_0..y_3,Degrees=>{{1,0},{1,0},{1,0},{0,1},{0,1},{0,1},{0,1}},Heft=>{1,1}]
   I = ideal(x_0^2+7*x_1^2+9*x_2^2,y_0*y_2*y_1^2,x_0*x_1)
   time Segre(I)
///
----------------
--Intro multi-proj Ex
--
----------------
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(4,CoefficientRing =>ZZ/32749,Variable =>symbol x)**projectiveSpace(2,CoefficientRing =>ZZ/32749,Variable =>symbol y)
    R=ring X
    gens R
    I=ideal(5*R_0*R_5,9*R_7*R_6*R_2-4*R_7^2*R_1)
    degrees I
    segre=time Segre(I)
    csm=time CSM(I)
    time Euler(I)
    time Euler(csm)

///
---------------------------------
--Multi-Proj Csm
-----------------------
--Ex. 1
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk);
    R=ring X;
    I=ideal(random({1,1},R),R_0^2*R_5^2-R_1*R_2*R_4*R_5)
    time CSM(I)

///

--Ex. 2
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(6,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk);
    R=ring X;
    I=ideal(R_0^2*R_1-R_2^3,R_7^2)
    degrees I
    codim I
    time CSM(I)

///

--Ex. 3
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(5,CoefficientRing =>kk)**projectiveSpace(3,CoefficientRing =>kk);
    R=ring X;
    I=ideal(4*R_0*R_6-7*R_7*R_2,R_0*R_4*R_8)
    degrees I
    codim I
    time CSM(I)

///

--Ex. 4
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(3,CoefficientRing =>kk);
    R=ring X;
    I=ideal((R_0*R_1-R_2^2)*R_4,R_5*(R_6^2-R_7*R_6))
    codim I
    time CSM(I)

///

--Ex. 5
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(3,CoefficientRing =>kk);
    R=ring X;
    I=ideal((R_0*R_1-R_2^2)*R_4,R_5*(R_6^2-R_7*R_6),R_0*R_3^2)
    K=ideal(I_0*I_1*I_2)
    codim I
    time CSM(I)
    time CSM(K)

///





------------------------------------------
--
--CSM complete intersection examples
-----------------------------------------

--Ex. 1
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk);
    R=ring X;
    A=ChowRing(R);
    I=ideal(random({1,1},R),random({1,1},R),R_1*R_0*R_3-R_0^2*R_4);
    csm=time CSM(A,I)
    csm2=time CSM(A,I,Method=>DirectCompleteInt)
    

///


--Ex. 2
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(3,CoefficientRing =>kk);
    R=ring X;
    A=ChowRing(R);
    I=ideal(random({2,1},R),R_1*R_0*R_4);
    csm=time CSM(A,I)
    csm2=time CSM(A,I,Method=>DirectCompleteInt)
///

--Ex. 3
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk);
    R=ring X
    A=ChowRing(R);
    I=ideal(random({2,1,0},R),random({0,0,1},R),R_2*R_6-7*R_0*R_7);
    csm=time CSM(A,I)
    csm2=time CSM(A,I,Method=>DirectCompleteInt)
    
   
///
--Ex. 4
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    kk=ZZ/32749;
    X=projectiveSpace(3,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk)**projectiveSpace(2,CoefficientRing =>kk);
    R=ring X;
    A=ChowRing(R);
    I=ideal(random({2,1,0},R),R_2*R_5-7*R_0*R_6);
    csm=time CSM(A,I)
    csm2=time CSM(A,I,Method=>DirectCompleteInt)
    
   
///
----------------------------------
--End csm examples
--
----------------------------------

----------------------------------
--Tests of multi-projective Segre
--
---------------------------------

--EX. 1
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(2,CoefficientRing =>ZZ/32749)**projectiveSpace(3,CoefficientRing =>ZZ/32749);
    R=ring X;
    I2 = ideal(random({2,0},R),R_3*R_5*R_4^2,R_0*R_1)
    time Segre(X,I2);
///

--EX. 2
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(1,CoefficientRing =>ZZ/32749)**projectiveSpace(1,CoefficientRing =>ZZ/32749)**projectiveSpace(1,CoefficientRing =>ZZ/32749);
    R=ring X;
    I2 = ideal(random({2,1,1},R),R_3*R_5*R_4^2-R_2*R_5^3)
    s=time Segre(X,I2);
///

--EX. 3
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(3,CoefficientRing =>ZZ/32749)**projectiveSpace(2,CoefficientRing =>ZZ/32749);
    S=ring X;
    I = ideal(S_0*S_1*S_2-S_2^2*S_3,S_0*S_2*S_1*S_3)
    time Segre(X,I)
///

--EX. 4
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(5,CoefficientRing =>ZZ/32749)**projectiveSpace(3,CoefficientRing =>ZZ/32749);
    S=ring X;
    I = ideal(S_0^2*S_7^2-S_6*S_7*S_1*S_2)
    degrees I
    time Segre(X,I)
///

--EX. 5
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(2,CoefficientRing =>ZZ/32749)**projectiveSpace(3,CoefficientRing =>ZZ/32749)**projectiveSpace(1,CoefficientRing =>ZZ/32749);
    S=ring X;
    I = ideal(S_0*S_1*S_3-5*S_4*S_2^2,S_4*S_6)
    degrees I
    time Segre(X,I)
///

--EX. 6
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(2,CoefficientRing =>ZZ/32749)**projectiveSpace(2,CoefficientRing =>ZZ/32749)**projectiveSpace(2,CoefficientRing =>ZZ/32749);
    S=ring X;
    I = ideal(S_0*S_3*S_6,S_5*S_7-7*S_4*S_8)
    degrees I
    time Segre(X,I)
///

--EX. 7
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(4,CoefficientRing =>ZZ/32749)**projectiveSpace(3,CoefficientRing =>ZZ/32749)**projectiveSpace(3,CoefficientRing =>ZZ/32749);
    S=ring X;
    gens S
    I = ideal(S_0*S_6^2+8*S_1*S_5*S_7,S_10*S_9*S_0-S_11^2*S_1)
    degrees I
    time Segre(X,I)
///

--EX. 8
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(4,CoefficientRing =>ZZ/32749)**projectiveSpace(3,CoefficientRing =>ZZ/32749)**projectiveSpace(5,CoefficientRing =>ZZ/32749);
    S=ring X;
    I = ideal(S_0*S_5^2-S_2*S_5*S_6,S_8*S_13-6*S_5*S_14)
    degrees I
    time Segre(X,I)
///

--EX. 9
TEST ///
{*
    restart
    needsPackage "MultiProjChar"
    needsPackage "NormalToricVarieties"
*}
    X=projectiveSpace(2,CoefficientRing =>ZZ/32749)**projectiveSpace(2,CoefficientRing =>ZZ/32749)**projectiveSpace(1,CoefficientRing =>ZZ/32749);
    S=ring X;
    gens S
    I = ideal(S_3*(S_4-S_5),S_5*S_4*S_6+9*S_3^2*S_7,S_0^2*S_4,S_0^2+S_1^2)
    degrees I
    time Segre(X,I)
///
--------------
--End Segre Exs
--
-------------