--*********************************************************************
--*
--*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
--*
--*********************************************************************

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

getrings := 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));

--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:=symbol t;
u:=symbol u;
z:=symbol z;
T:=kk[t_1 .. t_m];
U:=kk[t_1 .. t_m , z_0 .. z_n,MonomialOrder => Eliminate m];
irrel:=ideal(t_1 .. t_m);
irrelz:=ideal(z_0 .. z_n);
V:=kk[u,t_1 .. t_m,z_0 .. z_n,MonomialOrder => Eliminate 1];
M:=gens substitute(center,V);
N:=matrix{{t_1..t_m}}-u*M;
K:=gens gb ideal(N);
blowup:=ideal(selectInSubring(1,K));
blowup=substitute(blowup,U);
blowup=saturate(blowup,irrel);
d:=0;
Ide:=symbol Ide;Ide_0;jde:=symbol jde;
Ide#0=blowup;
while d<n do
 (d=d+1;Ide#d=saturate(Ide_(d-1)+substitute(ideal(random(1,T)),U),irrel); 
  if dim Ide_d-2 > n-d then (d=d-1, <<"."));
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 -> (
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;
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;
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:=symbol sset;sset_0;psum:=symbol psum;psum_0;
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;
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:=symbol sset;sset_0;psum:=symbol psum;psum_0;
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:=symbol sset;sset_0;psum:=symbol psum;psum_0;
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)
);

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