setRandomSeed 1
--*********************************************************************
--*
--*This file contains the code implementing the algorithms described
--*in my paper "Computing characteristic classes of projective schemes".
--*This code can be freely used, provided use is properly acknowledged.
--*Please send comments to    aluffi@math.fsu.edu
--*
--*Paolo Aluffi, April 2002
--*
--*********************************************************************

--*********************************************************************
--*
--*Minor changes implemented in January 2003:
--*extracted a separate routine blup for blow-ups, and replaced a dimension
--*check in presegre with a direct verification that an element is not
--*a zero-divisor.
--*
--*********************************************************************

--*********************************************************************
--*
--*Adapted to work with 0.9.95, June 2007
--*
--*********************************************************************

--*********************************************************************
--*
--*Adapted to work with 1.1, March 2008
--*
--*(avoiding an unnecessary transfer of a local variable, which was
--*confusing 1.1)
--*
--*********************************************************************

--getrings extracts the ring of the ideal and manufactures
--local copies of the coefficient field (kk), the ring (S, assumed 
--to be a polynomial ring), the dimension, and the ideal itself.

getrings := Icenter -> (
print Icenter;
S:=ring Icenter;
if not (isPolynomialRing S) then 
<<"Sorry, expecting an ideal defined over a polynomial ring"<<endl
else 
kk:=coefficientRing S;
n:=numgens S-1;
z:=symbol z;
S=kk[z_0..z_n];
sequence(kk,S,n,(map(S,ring Icenter,{z_0..z_n})) Icenter)
);

--fixideal prepares the ideal for use by (pre)segre. This amounts to
--trimming, fixing generators so that they all have the same 
--degree, and replacing the zero ideal by the ideal generated by 0

fixideal := (S,n,Icenter) -> (
cr:=0;en:=0;
tem:=ideal 0_S;
center:=trim substitute(Icenter,S);
m:=numgens center;
if m!=0 then (
ma:=first max degrees center;
ls:=first entries gens center;
for i from 0 to m-1 do
	(en=ls_i,cr=ma-first degree en;
	 for j from 0 to n do tem=tem+ideal(en*((entries vars S)_0_j)^cr));
	 tem=trim tem);
sequence(ma,tem));

--blup computes an ideal for the Rees algebra of the given ideal,
--for use in presegre.

blup = (I,t) -> (
	m := numgens I;
	S := ring I;
	kk := coefficientRing S;
	R := kk[gens S, t_0 .. t_(m-1)];
	II := substitute(I,R);
	J := ideal apply(0..(m-2), i -> apply((i+1)..(m-1), 
             j -> (II_i*t_j-II_j*t_i)));
        if m==1 then J=ideal 0_R;
	saturate(J,II_0)
);

--presegre applies the previous routines, then computes the shadow
--of the blow-up along the given ideal. This is encoded in the sequence
--of degrees of the images of suitable loci. The output includes
--the max degree of a generator of the ideal.

presegre = (kk,S,n,Icenter) -> (
init:=fixideal(S,n,Icenter);
center:=init_1;
m:=numgens center;
T:=kk[v_1 .. v_m];
U:=kk[v_1 .. v_m , gens S ,MonomialOrder => Eliminate m];
irrel:=ideal(v_1 .. v_m);
blowup:=substitute(blup(center,v),U);
blowup=saturate(blowup,irrel);
Ide:=new MutableList from {0..n};
d:=0;L:=0;Iq:=0;
Ide#0=blowup;
while d<n do
	(L=substitute(ideal(random(1,T)),U); 
	     Iq=(Ide#d:L);
	 if not(isSubset(Iq,Ide#d)) then <<"."
		else (d=d+1;Ide#d=saturate(Ide#(d-1)+L,irrel)); 
        );
tem:={};
for d from 0 to n do 
 (jde=substitute(selectInSubring(1,gens gb Ide#d),S);
  tem=append(tem,degree ideal(jde)));
sequence(init_0,tem));

--segre uses the sequence obtained in presegre, and computes the
--push-forward of the Segre class and of the Fulton class. These are 
--elements of the globally defined ring intringPn, generated by the 
--hyperplane class H. The global variables segreclass and fultonclass
--carry the two classes, available for further manipulations.
--segreclass is displayed.

segre = Icenter -> (
t:=symbol t;
init:=getrings(Icenter);
kk:=init_0;
S:=init_1;
n:=init_2;
center:=init_3;
init=presegre(kk,S,n,center);
ma:=init_0;tem:=init_1;
H=symbol H;
intringPn=ZZ[H]/(H^(n+1));
poly:=sum(0..n,s->tem#s*H^s*(1+ma*H)^(n-s));
segreclass=1 - poly * sum(0..n,i->binomial(n+i,i)*(-ma*H)^i);
<<"Segre class : "<<segreclass<<endl;
fultonclass=(1+H)^(n+1)*segreclass;
);

--CF does the same as segre, but displays fultonclass

CF = Icenter -> (
init:=getrings(Icenter);
kk:=init_0;
S:=init_1;
n:=init_2;
center:=init_3;
t:=symbol t;
init=presegre(kk,S,n,center);
ma:=init_0;tem:=init_1;
H=symbol H;
intringPn=ZZ[H]/(H^(n+1));
poly:=sum(0..n,s->tem#s*H^s*(1+ma*H)^(n-s));
segreclass=1 - poly * sum(0..n,i->binomial(n+i,i)*(-ma*H)^i);
fultonclass=(1+H)^(n+1)*segreclass;
<<"Fulton class : "<<fultonclass<<endl;
);

--csm computes the Chern-Schwartz-MacPherson class of a hypersurface,
--for use by the more general CSM. The input consists of the intersection
--ring W, where the answer is to be placed, the equation hyper of the
--hypersurface, and the information kk,S,n to be carried along to presegre.

csm := (kk,S,n,W,hyper) -> (
jac:=ideal jacobian ideal hyper;
t:= symbol t;
tem:=(presegre(kk,S,n,jac))_1;
--use W;
(1+H)^(n+1) - sum(0..n,d->tem#d*(-H)^d*(1+H)^(n-d))
);

--CSM assembles many csm computations to get the Chern-Schwartz-MacPherson
--class of an arbitrary ideal. The global variable csmclass carries the
--answer.

CSM = idea -> (
init:=getrings(idea);
kk:=init_0;
S:=init_1;
n:=init_2;
schem:=init_3;
H=symbol H;
intringPn=ZZ[H]/(H^(n+1));
r:=numgens schem;
sset:=new MutableList from {0..r};
psum:=new MutableList from {0..r};
gschem:=(entries gens schem)_0;
for s from 1 to r do
   (--<<";"<<endl;
   sset#s=apply(subsets(gschem,s),eq->csm(kk,S,n,intringPn,product(eq)));
   psum#s=sum(sset#s));
csmclass=sum(1..r,s->-(-1)^s*psum#s);
<<"Chern-Schwartz-MacPherson class : "<<csmclass<<endl;
);

--milnor does both CSM and CF, and defines a global variable milnorclass
--as well.

milnor = Icenter -> (
init:=getrings(Icenter);
kk:=init_0;
S:=init_1;
n:=init_2;
center:=init_3;
t:=symbol t;
init=presegre(kk,S,n,center);
ma:=init_0;tem:=init_1;
H=symbol H;
intringPn=ZZ[H]/(H^(n+1));
poly:=sum(0..n,s->tem#s*H^s*(1+ma*H)^(n-s));
segreclass=1 - poly * sum(0..n,i->binomial(n+i,i)*(-ma*H)^i);
fultonclass=(1+H)^(n+1)*segreclass;
<<"Fulton class : "<<fultonclass<<endl;
schem:=center;
r:=numgens schem;
sset:=new MutableList from {0..r};
psum:=new MutableList from {0..r};
gschem:=(entries gens schem)_0;
for s from 1 to r do
   (--<<";"<<endl;
   sset#s=apply(subsets(gschem,s),eq->csm(kk,S,n,intringPn,product(eq)));
   psum#s=sum(sset#s));
csmclass=sum(1..r,s->-(-1)^s*psum#s);
<<"Chern-Schwartz-MacPherson class : "<<csmclass<<endl;
milnorclass=csmclass-fultonclass;
<<"Milnor class : "<<milnorclass<<endl;
);

--euleraffinehyp computes the Euler characteristic of an affine 
--hypersurface, for use by euleraffine.

euleraffinehyp := (kk,S,n,W,eqn) -> (
enleqn:=homogenize(eqn,(entries vars S)_0_0)*(entries vars S)_0_0;
cuteqn:=substitute(enleqn,{(entries vars S)_0_0=>0});
tem:=csm(kk,S,n,W,enleqn)-csm(kk,S,n,W,cuteqn);
tem_(H^n)+1
);

--euleraffine computes the Euler characteristic of an affine scheme,
--by performing many euleraffinehyp.

euleraffine = idea -> (
S:=ring idea;
if not (isPolynomialRing S) then 
<<"Sorry, expecting an ideal defined over a polynomial ring"<<endl
else 
kk:=coefficientRing S;
n:=numgens S;
H=symbol H;
intringPn=ZZ[H]/(H^(n+1));
z:=symbol z;
S=kk[z_0..z_n];
schem:=(map(S,ring idea,{z_1..z_n})) idea;
r:=numgens schem;
gschem:=(entries gens schem)_0;
sset:=new MutableList from {0..r};
psum:=new MutableList from {0..r};
for s from 1 to r do
 (--<<";"<<endl;
 sset#s=apply(subsets(gschem,s),
             eq->euleraffinehyp(kk,S,n,intringPn,product(eq)));
 psum#s=sum(sset#s));
sum(1..r,s->-(-1)^s*psum#s)
);

--*
--*********************************************************************

--*********************************************************************
--*
--* Adding the following `test' code February 2009; it runs through
--* the examples in the documentation on 
--*
--* http://www.math.fsu.edu/~aluffi/CSM/CSMexamples.html
--*
--* and should raise an error message for any computation that does 
--* not produce the expected result.
--*
--*********************************************************************

-- (
       ringP2=QQ[x,y,z];
       ringP3=QQ[x,y,z,w];
       ringP4=QQ[x,y,z,w,t];
       use ringP3;
       CSM ideal 0_ringP3;
       assert(csmclass==4*H^3+6*H^2+4*H+1);
       threelines=ideal(x*y,x*z,y*z);
       segre threelines;
       assert(segreclass==-10*H^3+3*H^2);
       threecoplanarlines=ideal(z,x*y*(x+y));
       segre threecoplanarlines;
       assert(segreclass==-12*H^3+3*H^2);
       use ringP2; threecoplanarlines2=ideal(x*y*(x+y));
       segre threecoplanarlines2;
       assert(segreclass==-9*H^2+3*H);
       CF threelines;
       assert(fultonclass==2*H^3+3*H^2);
       CF threecoplanarlines;
       assert(fultonclass==3*H^2);
       CSM threelines;
       assert(csmclass==4*H^3+3*H^2);
       CSM threecoplanarlines;
       assert(csmclass==4*H^3+3*H^2);
       use ringP3;threeplanes=ideal(x*y*z);cubic=ideal(x^3+y^3+z^3+w^3);
       CF threeplanes;
       assert(fultonclass==9*H^3+3*H^2+3*H);
       CF cubic;
       assert(fultonclass==9*H^3+3*H^2+3*H);
       CSM threeplanes;
       assert(csmclass==4*H^3+6*H^2+3*H);
       CSM cubic;
       assert(csmclass==9*H^3+3*H^2+3*H);
       use QQ[x,y,z]; eu=euleraffine ideal(x^3+y^3+z^3-1);
       assert(eu==9);
       eu:=euleraffine ideal(x^3+y^3+z^2-1);
       assert(eu==5);
       use ringP3; twolinesnilp=ideal(x,z)*ideal(y,z);
       CF twolinesnilp; 
       assert(fultonclass==4*H^3+2*H^2);
       CSM twolinesnilp;
       assert(csmclass==3*H^3+2*H^2);
       twolinesred=ideal(x*y,z);
       CF twolinesred; 
       assert(fultonclass==2*H^3+2*H^2);
       CSM twolinesred;
       assert(csmclass==3*H^3+2*H^2);
       use ringP4; CSM ideal(x^5+y^5+z^5+w^5+t^5);
       assert(csmclass==-200*H^4+50*H^3+5*H);
       CSM ideal(x^3*t^2+x*z^4+w^5-y^2*t^3);
       assert(csmclass==4*H^4+38*H^3+5*H);
       use ZZ/2[x,y,z,w]; 
       q0 = 215*x^3*z^2+27*x^4*z+x^5+669*x^2*z^3+720*x*z^4;
       q1 = 5*x^4*y+27*x^4*w+108*x^3*y*z+430*x^3*z*w+645*x^2*y*z^2+720*y*z^4+1338*y*z^3*x+2007*x^2*z^2*w+2880*x*z^3*w;
       q2 = 10*x^3*y^2+215*x^3*w^2+108*x^3*y*w+162*x^2*y^2*z+1290*x^2*y*z*w+645*x*y^2*z^2+2007*x^2*w^2*z+669*y^2*z^3+2880*y*z^3*w+4014*x*y*z^2*w+4320*x*z^2*w^2;
       q3 = 10*x^2*y^3+215*y^3*z^2+162*x^2*y^2*w+645*x^2*y*w^2+108*x*y^3*z+1290*x*y^2*z*w+669*x^2*w^3+2880*x*w^3*z+2007*y^2*z^2*w+4014*x*y*z*w^2+4320*y*z^2*w^2;
       q4 = 5*x*y^4+27*y^4*z+108*x*y^3*w+645*x*y^2*w^2+430*y^3*z*w+720*x*w^4+1338*x*w^3*y+2007*y^2*w^2*z+2880*y*w^3*z;
       q5 = 215*y^3*w^2+27*y^4*w+y^5+669*y^2*w^3+720*y*w^4;
       dtupleideal= ideal(q0,q1,q2,q3,q4,q5);
       segre dtupleideal;
       assert(segreclass==-70*H^3+13*H^2);
       use ZZ/3[x,y,z,w]; 
       q0 = 215*x^3*z^2+27*x^4*z+x^5+669*x^2*z^3+720*x*z^4;
       q1 = 5*x^4*y+27*x^4*w+108*x^3*y*z+430*x^3*z*w+645*x^2*y*z^2+720*y*z^4+1338*y*z^3*x+2007*x^2*z^2*w+2880*x*z^3*w;
       q2 = 10*x^3*y^2+215*x^3*w^2+108*x^3*y*w+162*x^2*y^2*z+1290*x^2*y*z*w+645*x*y^2*z^2+2007*x^2*w^2*z+669*y^2*z^3+2880*y*z^3*w+4014*x*y*z^2*w+4320*x*z^2*w^2;
       q3 = 10*x^2*y^3+215*y^3*z^2+162*x^2*y^2*w+645*x^2*y*w^2+108*x*y^3*z+1290*x*y^2*z*w+669*x^2*w^3+2880*x*w^3*z+2007*y^2*z^2*w+4014*x*y*z*w^2+4320*y*z^2*w^2;
       q4 = 5*x*y^4+27*y^4*z+108*x*y^3*w+645*x*y^2*w^2+430*y^3*z*w+720*x*w^4+1338*x*w^3*y+2007*y^2*w^2*z+2880*y*w^3*z;
       q5 = 215*y^3*w^2+27*y^4*w+y^5+669*y^2*w^3+720*y*w^4;
       dtupleideal= ideal(q0,q1,q2,q3,q4,q5);
       segre dtupleideal;
       assert(segreclass==-58*H^3+11*H^2);
       use ringP2; milnor ideal(y^6+z*x^3*y^2+z^2*x^4);
       assert(fultonclass==-18*H^2+6*H);
       assert(csmclass==6*H);
       assert(milnorclass==18*H^2);
       use ringP3; milnor ideal(y^2*z-x^3-x^2*z);
       assert(fultonclass==9*H^3+3*H^2+3*H);
       assert(csmclass==2*H^3+4*H^2+3*H);
       assert(milnorclass==-7*H^3+H^2);
-- )
   
--*********************************************************************