words(0):kernelopts(printbytes=false):

# At the bottom: an attempt to program Rahul's relations of 1999-01-12.

lam:=proc(n)
# much easier!
option remember:
expand(coeff(taylor(exp(add(bernoulli(2*i)*ka||(2*i-1)*t^(2*i-1)/2/i/(2*i-1),
i=1..iquo(n+1,2))),t,n+1),t,n)):
end:

lamb:=proc(n)
# as above, with ch classes
# Much easier, so I replace "lambda" with "lamb" in rest of file.
# Only occurrence: in "klp".
option remember:
expand(coeff(taylor(exp(add(ch||(2*i-1)*t^(2*i-1)*(2*i-2)!,
i=1..iquo(n+1,2))),t,n+1),t,n)):
end:

oddlist:=proc(n)
# superfluous now, see below.
local k:
if n<1 then [] else
k:= iquo(n+1,2):
hulp(k,n)
fi:
end:

hulp:=proc(k,n)
##Generates a list of [d1,...,dk] such that 0<=di<=n,
## sum((2*i-1)*di,i=1..k)=n
# superfluous now, see below.
local a,i,j,result:
if k=1 then
        if n<0 then [] else [[n]] fi:
else
        result:=[]:
        for i from n by -1 to 0 do
        a:=hulp(k-1,n-(2*k-1)*i):
        result:=[op(result),[op(a[j]),i]$j=1..nops(a)]
        od:
result;
fi:
end:
 
lambda:=proc(n)
# superfluous now, "lamb" is easier.
local class,odd,i,elt,term,j,a;
class:=0:
odd:=oddlist(n):
if n=0 then class:=1 fi:
for i to nops(odd) do
elt:=op(i,odd):
term:=1:
for j to nops(elt) do
a:=op(j,elt):
if a>0 then term:=term*((2*j-2)!*ch||(2*j-1))^a/(a!) fi:
od:
class:=class+term:
od:
expand(class):
end:

read `KT21v5.m`:
read `tauZaalV5.m`:	## Zaal's work.
read `minderV5.m`:

klmgn:=proc(g,n,y)  	## g is the genus; y is a list of numbers of psi's,
  			## kappa's and ch_i's on \bar{M}_{g,n} . The procedure
   			## returns the intersection number.
option remember:
local aa,x,nx,ex,ex2,sx,ans,pp,npp,b,i2,bb,m,facm,h,pol,npol,temp,bh,bgh:
aa:=4*g-3+2*n:  ## Number of classes considered (psi's, kappa's, odd ch_i's).
x:=sort(y):
nx:=nops(x):
ex:=expo(x):
ex2:=[op(ex),0$(aa-nops(ex))]:
sx:=add(ex2[i],i=1..n)
   +add((i-n)*ex2[i],i=n+1..3*g-3+2*n)
   +add((2*(i-3*g+3-2*n)-1)*ex2[i],i=3*g-2+2*n..aa):
## Weighted sum of entries should be 3g-3+n.
if not sx=3*g-3+n then [g,n,y] elif nx=0 then 1 else
 if x[nx]<3*g-2+2*n 		# if 1 (only psi's, kappa's)
 then ans:=taug(g,n,x)  	# This includes g=0.
 else 
  if n>0 and g>1
  then pp:=product('(ps||i)^ex2[i]','i'=1..n) *
           product('(ka||(i-n))^ex2[i]','i'=n+1..3*g-3+2*n):
  pp:=push(g,n,pp):
   if type(pp,`+`)
   then npp:=nops(pp): for i2 from 1 to npp do b[i2]:=op(i2,pp) od:
   else npp:=1: b[1]:=pp:
   fi:
  ans:=0:
if not b[1]=0 then			### new (2 nov. 1998)
   for i2 from 1 to npp do
   bb:=['degree(b[i2],ka||j)'$'j'=1..3*g-3]:
   ans:=ans+
   lcoeff(b[i2])*klmgn(g,0,[''j'$bb[j]'$'j'=1..3*g-3,
   ''j-2*n'$ex2[j]'$'j'=3*g-2+2*n..aa]):
   od:
fi:					### new (2 nov. 1998)
  elif g=1 then   		# Get rid of that single lambda there.
if ex2[2*n+1]=1 then
  ans:=1/24*klmgn(0,n+2,
[''j'$ex2[j]'$'j'=1..n,''j+2'$ex2[j]'$'j'=n+1..2*n]):
else ans:=0 fi:
else     		# we are now on Mgbar with g at least 2, with a ch_i.
# find the last index, rewrite that ch_i according to Mumford's formula||
# That is pretty much it. Except that I probably want to push-down to Mhbar.
m:=2*(x[nx]-3*g+3)-1:  		## the component of the Chern character.
facm:=bernoulli(m+1)/(m+1)!: 	## the factor.
#ans:=klmgn(g,0,['x[i]'$'i'=1..nx-1,m]):  	## the kappa term.
ans:=klmgn(g,0,sort(['x[i]'$'i'=1..nx-1,m])):
# now the Delta_0 term.
# should go only to m here, and take one m away.
ans:=ans+
1/2*add((-1)^j*
klmgn(g-1,2,[1$j,2$(m-1-j),''i+2'$ex2[i]'$'i'=1..3*g-4,
''i+1'$ex2[i]'$'i'=3*g-2..(x[nx]-1),(x[nx]+1)$(ex2[x[nx]]-1)]),j=0..m-1):
# now the reducible terms.
for h from 1 to g-1 do
# write out the relevant terms, then apply klmgn
pol:=add((-D1)^j*F1^(m-1-j),j=0..m-1):
pol:=pol*subs(['D||i=0'$'i'=3*h..3*g-2,'F||i=0'$'i'=3*(g-h)..3*g-2],
product('(D||(i+1)+F||(i+1))^ex2[i]','i'=1..3*g-3)):
pol:=pol*subs(['D||i=0'$'i'=4*h..aa,'F||i=0'$'i'=4*(g-h)..aa],
# go only to m, and take m with one fewer.
product('(D||(i-3*g+2+3*h)+F||(i-3*g+2+3*(g-h)))^ex2[i]','i'=3*g-2..x[nx]-1)
*(D||(x[nx]-3*g+2+3*h)+F||(x[nx]-3*g+2+3*(g-h)))^(ex2[x[nx]]-1)):
pol:=expand(pol):
if type(pol,`+`) then npol:=nops(pol): for i2 from 1 to npol
do b[i2]:=op(i2,pol) od: else npol:=1: b[1]:=pol: fi:
temp:=0:
for i2 from 1 to npol do
bh:=['degree(b[i2],D||i)'$'i'=1..4*h-1]:
bgh:=['degree(b[i2],F||i)'$'i'=1..4*(g-h)-1]:
if bh[1]+add(bh[i]*(i-1),i=2..3*h-1)
+add(bh[i]*(2*(i-3*h+1)-1),i=3*h..4*h-1) = 3*h-2
then temp:=temp+lcoeff(b[i2])*
klmgn(h,1,['i$bh[i]'$'i'=1..4*h-1])*
klmgn(g-h,1,['i$bgh[i]'$'i'=1..4*(g-h)-1]):
fi:
od:
ans:=ans+1/2*temp:
od:    			## end of "for h from 1 to g-1 do"
ans:=ans*facm:
  fi:   							# fi 1
 fi:								# fi 0
#[ans,factor(pol)]:	## alternative for trouble-shooting purposes.
ans:
fi:
end:

taug:=proc(g,n,x)	## g is the genus; x is an ordered list of numbers
    			## of divisors and classes on \bar{M}_{g,n}, of weight
   			## 3*g-3+n. The numbers are all < 3g-2+2n.
   			## This corresponds to \psi's and \kappa's.
   			## The procedure returns the intersection number.
option remember:
local nx,ex,ex1,ex2,ans,j,i,kap,tu,np,a,pp,spp,part,j2,j3:
nx:=nops(x):
ex:=expo(x):
ex1:=[0$(n-nops(ex)),op(ex)]:
ex2:=sort(ex1):
if nx=0 or x[nx]<n+1 then       ## Only \psi's.
ans:=tau(ex2): 
else    			## Also \kappa's.
pp:=[''i'$ex[n+i]'$'i'=1..nops(ex)-n]:
spp:=add(pp[i],i=1..nops(pp)):
part:=pa[spp]:			### spp < 26; read `minder.m`.
for j from 1 to nnp[spp] do
if part[j]=pp then j2:=j: j:=nnp[spp] fi od:
tu:=KT[spp+2][j2]:		## spp <= 19 for now.
ans:=0:
if type(tu,`+`) then
np:=nops(tu):
for j from 1 to np do a[j]:=op(j,tu) od:
for j from 1 to np do ans:=ans+lcoeff(a[j])*
tau(sort(['ex1[i]'$'i'=1..n,''i'$degree(a[j],ta||i)'$'i'=2..3*g-2+n])):
od:
else ans:=lcoeff(tu)*
tau(sort(['ex1[i]'$'i'=1..n,''i'$degree(tu,ta||i)'$'i'=2..3*g-2+n])):
fi:
fi:
ans:
end:

## the procedure expo counts how many ones, twos, etc., the list x
## of positive integers contains.
expo:=proc(x)
local res,lijst,lang,big,j,i,i2:
if x=[] then [] else
lijst:=sort(x):
lang:=nops(lijst):
big:=lijst[lang]:
for j from 1 to big do res[j]:=0 od:
for i from 1 to lang do
res[x[i]]:=res[x[i]]+1:
od:
['res[i2]'$'i2'=1..big]:
fi:
end:

tau:=proc(x)
## x is the list [d_1, ..., d_n] ; d_i <= d_j als i <= j .
option remember:
local n,g,i,j,ans:
n:=nops(x):
g:=add(x[i]-1,i=1..n)/3+1:
if g=0 then ans:=(n-3)!/product('x[i]!','i'=1..n)
elif x[1]=0 then ans:=0: for j from 2 to n do if x[j]>0 then
ans:=ans+tau(sort(['x[i]'$'i'=2..j-1,x[j]-1,'x[i]'$'i'=j+1..n])) fi od:
elif x[1]=1 then if g=1 and n=1 then ans:=1/24 else
ans:=(2*g-2+n-1)*tau(['x[i]'$'i'=2..n]) fi:
else if type(g,integer) then ans:=TAU[x]: else ans:=0 fi:
fi:
ans:
end:

klrij:=proc(g,n)  	## list of relevant classes on \bar{M}_{g,n} .
option remember:
['ps||i'$'i'=1..n,'ka||i'$'i'=1..3*g-3+n,'ch||(2*i-1)'$'i'=1..g]:
end:

pushn:=proc(g,n,p)
local d,q,j,q0,nq,k,a:
d:=3*g-3+n:
q:=p:
for j from 1 to d do
q:=subs(ka||j=ka||j+(ps||n)^j,q)
# writing kappa's on level n in terms of kappa's pulled-back from level n-1.
od:
q0:=expand(subs(ps||n=0,q)):
q:=q-q0:
q:=collect(q,ps||n):
for j from d by -1 to 2 do
q:=subs((ps||n)^j=ka||(j-1),q)
od:
q:=subs(ps||n=2*g-3+n,q):
if type(q0,`+`) then nq:=nops(q0):
for k from 1 to nq do a[k]:=op(k,q0) od:
else nq:=1: a[1]:=q0:
fi:
for k from 1 to nq do
for j from 1 to n-1 do
if subs(ps||j=0,a[k])=0 then q:=q+expand(a[k]/ps||j) fi
# string equation.
od od:
q:
end:

push:=proc(g,n,p)
local q,j:
q:=p:
for j from n by -1 to 1 do q:=pushn(g,j,q) od:
expand(q):
end:

sub:=proc(g,n,p)
local dp,i,np,j,a,test,pp:
if type(p,`+`) then np:=nops(p): for j from 1 to np do a[j]:=op(j,p) od:
else np:=1: a[1]:=p: fi:
test:=0:
for j from 1 to np do
pp:=a[j]:
dp:=add(degree(pp,D||i),i=1..n)
   +add((i-n)*degree(pp,D||i),i=n+1..3*g-3+2*n)
   +add((2*(i-3*g+3-2*n)-1)*degree(pp,D||i),i=3*g-2+2*n..4*g-3+2*n):
## Weighted sum of entries should be 3g-3+n.
if not dp=3*g-3+n then test:=[g,n,pp] else
test:=test+
lcoeff(a[j],['D||i'$'i'=1..4*g-3+2*n])*
klmgn(g,n,[''i'$degree(a[j],D||i)'$'i'=1..4*g-3+2*n]):
#lcoeff(a[j])*op(1,mgn(g,[''i'$degree(a[j],D||i)'$'i'=1..(2^(n-1)*(g+1)+1)])):
#       alternative for trouble-shooting purposes.
fi:
od:
test:
end:

klp:=proc(g,n,x)
local ef:
ef:=subs(['la||i=lamb(i)'$'i'=1..g],x):
ef:=subs(['ch||(2*i-1)=D||(3*g-3+2*n+i)'$'i'=1..g],ef):
ef:=subs(['ps||i=D||i'$'i'=1..n],ef):
ef:=subs(['ka||i=D||(i+n)'$'i'=1..3*g-3+n],ef):
ef:=expand(ef):
sub(g,n,ef):
end:

# An attempt to program Rahul's relations of 1999-01-12.

rel:=proc(g,d,a)
# g>=2, d>=1, a=[a1,...,a_n] n non-negative numbers, with n>=1.
# Outcome: a relation in the tautological ring of M_g in dimension
# g-2+d - \sum a_i , so in codimension 2g-1 -d +\sum a_i . So
# we need d - \sum a_i >= g+1 to get possibly non-vanishing relations;
# but relations in codim.>= g-1 are interesting also, as they could
# provide an alternative proof of Looijenga's vanishing.
local pd,rel,i,m,lm,em,autm,facm,n,codim,x:
pd:=pa[d]:
rel:=0:
n:=nops(a):
for i from 1 to nnp[d] do
m:=pd[i]:
lm:=nops(m):
em:=expo(m):
autm:=product('em[z]!','z'=1..nops(em)):
facm:=(-1)^lm*product('m[z]^(m[z]-1)/m[z]!','z'=1..lm)/autm:
codim:=2*g-1-d+n+lm:
if min(0,codim)=0 then
x:=coeff(taylor(product('1/(1-m[z]*ps||(n+z)*t)','z'=1..lm),t,codim+1),t,codim):
x:=x*product('(ps||z)^a[z]','z'=1..n):
rel:=rel+facm*push2(g,n+lm,x): 
fi:
od:
expand(rel):
end:

# I want an alternative to "push" since the relations are more naturally
# interpreted in tau-expressions. In any case this is so in the top degree
# 2g-3+n  i.e. dimension g. It's easier if we try to verify the conjecture
# for lambda_g times psi's.
# At the moment I'm a little confused about the meaning in higher dimension.

# a variant of pushn: only valid for monomials in psi's; get rid of tau0
# and tau1; use the tau-expressions as a basis downstairs instead of
# kappa-monomials. E.g. for degree 3 I would want something 
# like tau4, tau23, tau222.
 
push2:=proc(g,n,p)
local q,j,ans,nq,k,a,mon,pp,dm,spp,part,j2:
q:=p:
for j from n by -1 to 1 do q:=pushn(g,j,q) od:
q:=expand(q):
# this is now a polynomial in kappa's. Now we will undo the kappa's in
# tau's using the list KT. Hence degree 19 at most in kappa's
# (g<=11 for Mtilde).
ans:=0:
if type(q,`+`) then nq:=nops(q):
for k from 1 to nq do a[k]:=op(k,q) od:
else nq:=1: a[1]:=q:
fi:
for k from 1 to nq do
mon:=a[k]: 	# monomial in kappa's
pp:=[]:
for j from 1 to 2*g-3 do	   # For Mtilde_g .
dm:=degree(mon,ka||j):
if dm>0 then pp:=[op(pp),j$dm] fi:
od:
# pp is the monomial as a list of nonzero indices.
spp:=add(pp[i],i=1..nops(pp)): # should be independent of k.
part:=pa[spp]:                     # spp < 26; read `minder.m`.
for j from 1 to nnp[spp] do
if part[j]=pp then j2:=j: j:=nnp[spp] fi od:
ans:=ans+lcoeff(mon)*KT[spp+2][j2]:            # spp <= 19 for now.
od:
rew(g,ans):
end:

rew:=proc(g,q)
# ta2^2 -> ta[[2, 2]], etc.
local ans,nq,k,a,mon,pp,j,dm:
ans:=0:
if type(q,`+`) then nq:=nops(q):
for k from 1 to nq do a[k]:=op(k,q) od:
else nq:=1: a[1]:=q:
fi:
for k from 1 to nq do
mon:=a[k]:      # monomial in tau's.
pp:=[]:
for j from 2 to 2*g-2 do           # For Mtilde_g .
dm:=degree(mon,ta||j):
if dm>0 then pp:=[op(pp),j$dm] fi:
od:
# pp is the monomial as a list of nonzero indices.
ans:=ans+lcoeff(mon)*ta[pp]:
od:
ans:
end:

temp:=proc(g,a) klp(g,nops(a),la||g*product('ps||i^(a[i]+1)','i'=1..nops(a)))/
klp(g,1,la||g*ps1^(2*g-2)) end:

effe:=proc(a) (add(a[i]+1,i=1..nops(a)))!/
product('(a[i]+1)!','i'=1..nops(a)) end:

vv:=proc(g,y) map(x->x/lcoeff(x),{op(y)}) minus {ta[[2*g-2]]} end:

