# 	 MAPLE script for computing numbers such as  of 
#	        CONICS IN P3 INCIDENT TO A  LINE
#		AND CONTAINED IN A CUBIC SURFACE
# adapted by A. Meireles from a script by P. Meurer ,
# #The number of rational quartics on Calabi-Yau hypersurfaces in
# weighted
# #projective space $\ps{2,1^4}$} Math. Scandin, #78, 63-83 1996,
# #http://arXiv.org/abs/math/alg-geom/9409001
 n:=3;
 xx:= seq(x.i, i=0..n); 
 ww:= seq(w.i, i=0..n);
 sf:=proc(vars,d) local i;
   coeff(expand(series(product( (1/(1-t*vars[i])),
         i=1..nops(vars) ),t=0,d+1)),t^d)
 end:
 F:=sf([xx],1); S2F:=sf([xx],2); S3F:=sf([xx],3);
 # prodwts(H) product of all weights
 prodwts:= proc(H) local t, cof, mon, res, i; 
     cof:= [coeffs(H, [seq(x.i, i=0..n)],'mon' )]:
     mon:= subs(seq(x.i=t^w.i, i=0..n),[mon]):
     res:=1:
     for i from 1 to nops(mon) do
          res:= res*subs(t=1,diff(mon[i],t))^cof[i]:
 od: end:
 # sumwts(H) sum of weights
  sumwts:= proc(H) local t, cof, mon, res, i; 
     cof:= [coeffs(H, [seq(x.i, i=0..n)],'mon' )]:
     mon:= subs(seq(x.i=t^w.i, i=0..n), [mon]):
     res:=0:
     for i from 1 to nops(mon) do
         res:= res+subs(t=1,diff(mon[i],t))*cof[i]:
 od: end:
 # cont counts the number of fixpts
 cont:=0:
 for i  to nops(F) do
   h:=op(i,F):
   S2F_h:=expand(S2F-h*F):
   for j to nops(S2F_h) do
     q:=op(j,S2F_h):
     cont:=cont+1:
     T[cont]:=expand( (F-h)/h + (S2F_h-q)/q ):
     K[cont]:=1/h:
     Q[cont]:=1/q: 
     expand(h*S2F):
     expand(q*F):
     idd:={op("),op("")}:
     V[cont]:=expand( S3F - convert(idd,`+`) ):
 od: od:
 WG:=[-2,3,41,119]; 
 wg:= seq(WG[i], i=1..n+1);
 #  N(a,b,c,d)=deg(c1(K)^a.c1(Q)^b.(c1(Q)+2c1(K))^c.c7(V)^d),
 #    for a+b+c+7d=8
 N:=proc(a1,a2,a3,a4) local i,r,den,f1,f2,g,f;
  r:=0: 
  for i from 1 to cont do  
    den[i]:=unapply(prodwts(T[i]), ww):    
    den[i]:=den[i](wg):
    f1[i]:=unapply(sumwts(K[i]), ww):
    f1[i]:=f1[i](wg):
    f2[i]:=unapply(sumwts(Q[i]), ww):
    f2[i]:=f2[i](wg): 
    g[i]:=unapply(prodwts(V[i]), ww):
    g[i]:=g[i](wg):
    f[i]:=f2[i]+2*f1[i]: 
    r:=r+(f1[i]^a1*f2[i]^a2*f[i]^a3*g[i]^a4)/den[i]:
 od: end:
 print('N(3,5,0,0)'=N(3,5,0,0)): print('N(0,0,8,0)'=N(0,0,8,0)):
 print('N(1,0,7,0)'=N(1,0,7,0)): print('N(2,0,6,0)'=N(2,0,6,0)):
 print('N(3,0,5,0)'=N(3,0,5,0)): print('N(1,0,0,1)'=N(1,0,0,1)):
 print('N(0,0,1,1)'=N(0,0,1,1)):