# This program is by Lothar Gottsche
# Warning:  The degrees calculated are not that of just irreducible curves

with(share): readshare(SF,combinat): with(SF):
with(linalg);



partr:=proc(bb);#transforms partition into short notation
l:=nops(bb);erg:=[];
if l>0 then 
dd:=bb[1];
for n from 1 to dd do
a[n]:=0;od;
for i from 1 to l do
a[bb[i]]:=a[bb[i]]+1;
od;
erg:=[];
for i from 1 to dd do
erg:=[op(erg),a[i]];
od;
fi;
erg;
end;


#The routine below calculates the degree of the locus of  nodal curves with
#"del" nodes in the linear series of curves of degree d on P^2. 
#To use it, type U(d,del);

#If you want more detailed info, as in the Caparoso-Harris paper, you can
#get that too by typing N(d,del,[alpha],[beta]);  where alpha and beta are
#the lists of alphas and betas as in the paper.


#N^d,delta on P2
U:=proc(d,del);
N(d,del,[],[d]);
end;

#starting value for N
M:=proc(del,a,b);
MM:=0;
if del=0 and ((a=[] and b=[1]) or (a=[1]  and b=[])) then 
MM:=1 fi;
MM;
end;

 

#N^d[alpha,beta] on P_2
N:=proc(d,del,a,b) option remember;
if d=1 then NN:=M(del,a,b);
else
aa:= a; bb:= b;
na:=nops(aa); nb:=nops(bb);
uu:=na-nb;
for i from 1 to -uu do 
aa:=[op(aa),0]; od;
for i from 1 to uu do 
bb:=[op(bb),0]; od;
n:=nops(aa);


NN:=0;
for k from 1 to nb do
if bb[k]>0 then 
NN:=NN+k*N(d,del,compress(subsop(k=aa[k]+1,aa)),compress(subsop(k=bb[k]-1,bb)));
fi;
od;
for da from 0 to d-1 do
db:=d-1 -da;
for  xx in   SF[Par](da) do  
aa0:=partr(xx);
  if bigger(aa,aa0) then
for  yy in   SF[Par](db) do
bb0:=partr(yy);
 if bigger(bb0,bb) then
deldel:=del+betr(mins(bb0,bb))-d+1;
if deldel<=del and deldel>=0 then
NN:=NN+IH(mins(bb0,bb))*binom(aa,aa0)*binom(bb0,bb)*N(d-1,deldel,
compress(aa0),compress(bb0));
fi;fi;od;fi;od;od;
fi;
NN;
end;
  

#The routine below calculates degrees of the locus of curves of arithmetic
#genus g in the linear series (a,b).  To use it type UP(a,b,g); where 
#a, b, and g are as in the above description.  I'm not sure what the next
#one does, since I haven't see Ravi's algorithm.

#N^D,delta for  P1 x P1

UP:=proc(df,dg,g);
NP(df,dg,g,[],[df]);
end;

#Starting value for NP
MP:=proc(df,g,a,b);
MM:=0;
if df>0 and g=1-df then
for i from 0 to df do 
if a=compress([i]) and b=compress([df-i]) then 
MM:=1 fi;od;
fi;
#lprint([df,g,a,b],MM);
MM;
end;

#N^D,delta[alpha,beta] for P1 X P1
NP:=proc(df,dg,g,a,b) option remember;
d:=df;
if dg=0 then NN:=MP(df,g,compress(a),compress(b));
else
aa:= a; bb:= b;
na:=nops(aa); nb:=nops(bb);
uu:=na-nb;
for i from 1 to -uu do 
aa:=[op(aa),0]; od;
for i from 1 to uu do 
bb:=[op(bb),0]; od;
n:=nops(aa);


NN:=0;
for k from 1 to nb do
if bb[k]>0 then 
NN:=NN+k*NP(df,dg,g,compress(subsop(k=aa[k]+1,aa)),compress(subsop(k=bb[k]-1,bb)));
fi;
od;
for da from 0 to d do
db:=d -da;
for  xx in   SF[Par](da) do  
aa0:=partr(xx);
  if bigger(aa,aa0) then
for  yy in   SF[Par](db) do
bb0:=partr(yy);
 if bigger(bb0,bb) then
gg:=g-betr(mins(bb0,bb))+1;
 NN:=NN+IH(mins(bb0,bb))*binom(aa,aa0)*binom(bb0,bb)*NP(df, dg-1,gg,
compress(aa0),compress(bb0))
;fi;od;fi;od;od;
fi;
NN;
end;




#N^D,delta for rational ruled surface

#The routines below calculate degrees for ruled surfaces
#For the surface F_e, the picard group is spanned by the divisors F and G,
#with F the class of a fibre and G the unique curve with self-intersection -e.
#A linear system on F_e is specified as (df)F + (dg)G with df and dg
#integers.  The routine calculates the degree of the locus in the linear
#series of curves of arithmetic genus g.  The function is called with
#UR(e,df,dg,g);  

#The more sophisticated numbers from the recursion is availible from the
#procedure NR, called with NR(e,df,dg,g,[alpha],[beta]); where alpha and
#beta are lists of moving and fixed contact with the curve G.

UR:=proc(e,df,dg,g);
d:=df-e*dg;
NR(e,df,dg,g,[],[d]);
end;



#Starting value for NR

MR:=proc(df,g,a,b);
MM:=0;
if df>0 and g=1-df then
for i from 0 to df do 
if a=compress([i]) and b=compress([df-i]) then 
MM:=1 fi;od;
fi;
#lprint([df,g,a,b],MM);
MM;
end;

#N^D,delta[alpha,beta] for rational ruled surface
NR:=proc(e,df,dg,g,a,b) option remember;
#lprint(e,df,dg,g,a,b);
d:=df-e*dg;
if dg=0 then NN:=MR(df,g,compress(a),compress(b));
else
aa:= a; bb:= b;
na:=nops(aa); nb:=nops(bb);
uu:=na-nb;
for i from 1 to -uu do 
aa:=[op(aa),0]; od;
for i from 1 to uu do 
bb:=[op(bb),0]; od;
n:=nops(aa);


NN:=0;
for k from 1 to nb do
#lprint(k);
if bb[k]>0 then 
NN:=NN+k*NR(e,df,dg,g,compress(subsop(k=aa[k]+1,aa)),compress(subsop(k=bb[k]-1,bb)));
fi;
od;
for da from 0 to d+e do
db:=d+e -da;
#lprint([da,db]);
for  xx in   SF[Par](da) do  
aa0:=partr(xx);
  if bigger(aa,aa0) then
for  yy in   SF[Par](db) do
bb0:=partr(yy);
 if bigger(bb0,bb) then
gg:=g-betr(mins(bb0,bb))+1;
 NN:=NN+IH(mins(bb0,bb))*binom(aa,aa0)*binom(bb0,bb)*NR(e,df, dg-1,gg,
compress(aa0),compress(bb0))
;fi;od;fi;od;od;
fi;
NN;
end;




compress:=proc(a);
aa:=a;
n:=nops(aa); if n>0 then
if aa[n]=0 then aa:=compress([aa[1..n-1]]);fi; 
# else aa:=[0]; 
fi;
aa;
end;

binom:=proc(a,b);
aa:=a;bb:=b;
na:=nops(a);
nb:=nops(b);
uu:=na-nb;
for i from 1 to -uu do 
aa:=[op(aa),0]; od;
for i from 1 to uu do 
bb:=[op(bb),0]; od;
n:=nops(aa);
ee:=1;
for i from 1 to n do
ee:=ee*binomial(aa[i],bb[i]);
od;
ee;
end;



IH:=proc(a,b);
n:=nops(a);
ee:=1;
for i from 1 to n do
ee:=ee*i^a[i] od;
ee;
end;

II:=proc(a);
n:=nops(a);
ee:=0;
for i from 1 to n do ee:=ee+i*a[i]; od;
ee; end;

#|partition|
betr:=proc(a);
n:=nops(a);
ee:=0;
for i from 1 to n do ee:=ee+ a[i]; od;
ee; end;

#difference of partitions
mins:=proc(a,b);
aa:=a;bb:=b;
na:=nops(a);
nb:=nops(b);
uu:=na-nb;
for i from 1 to -uu do 
aa:=[op(aa),0]; od;
for i from 1 to uu do 
bb:=[op(bb),0]; od;
n:=nops(aa);cc:=[];
for i from 1 to n do
cc:=[op(cc),aa[i]-bb[i]]
od;
cc;
end;

#compare partitions
bigger:=proc(a,b);
cc:=mins(a,b);
n:=nops(cc);
ee:=1;
for i from 1 to n do
if cc[i]<0 then ee:=ee-1;fi;
od;
ee>0;
end;



#List of N^D,delta for P2
ww:=proc(n);
for d from 1 to n do
for del from 0 to binomial(d,2)+1 do
lprint([d,del], U(d,del));
od;od;
end;

#List of N^D,delta for P1 X P1
wwP:=proc(n,m);
for df from 1 to n do
for dg from 1 to m do 

for del from 0  to n*m+1 do
g:=n*m-n-m+1-del;
lprint([df,dg,del,g], UP(df,dg,g));
od;od;od;
end;



#list of N^d,delta[alpha,beta] for P2
vv:=proc(n);
for d from 1 to n do


for del from 0 to binomial(d,2)+1 do
for da from 0 to d do
db:=d-da;
 for  xx in   SF[Par](da) do  
aa:=partr(xx);
 for  yy in   SF[Par](db) do  
bb:=partr(yy);

lprint([d,del,aa,bb], N(d,del,aa,bb));
od;od;od;od;od;
end;


#Power series for P1 X P1
nPser:=proc(d1,d2, del);
AA:=array(1..d1,1..d2);
for n from 1 to d1 do
for m from n to d2 do

AA[n,m]:=0;
for dell from 0 to min(del,n*m+2) do
g:=n*m-n-m+1-dell;
AA[n,m]:=AA[n,m]+UP(n,m,g)*q^dell;
od;
lprint([n,m],series(AA[n,m],q,min(del+1,n*m+2)));
od;od;
end;


#Power series for ruled surface

Rser:=proc(e1,d1, del);
 
for e from 1 to e1 do
for n from 1 to d1 do
a:=round(n-e);

for m from 1 to a do

r:=n*m +n+m-e*binomial(m+1,2);
 if r>=0 then
AA :=0;
for dell from 0 to del do
g:=n*m-n-m-binomial(m,2)*e+1-dell;
 


AA :=AA +UR(e,n,m,g)*q^dell;

od;
lprint([e,n,m],series(AA,q,del+1));fi;
od;od;od;
end;