gc(0):
# Start by defining a few procedures that are convenient when 
# computing in the representation ring of a torus.

# setuptorus(n,m) sets up an (n+1)-dimensional torus with
# characters x0,...,xn. The representation ring is identified
# with the Laurent series ring Z[xi,1/xi, i=0..n]. The canonical
# representation sum(x.i,i=0..N) is denoted by V[1]. Its symmetric
# powers are denoted by V[k].

setuptorus:=proc(N,M) local genf,i,k;
	genf:=expand(convert(series(
		1/product('1-x.i*t',i=0..N),t,M+1),polynom));
	for i from 0 to M do V[i]:=sort(coeff(genf,t,i)) od; i:='i';
	dualsubstring:=seq(x.i=x.i^(-1),i=0..N);
	variables    :=[seq(x.i,i=0..N),seq(x.i^(-1),i=0..N)];
	ranksubstring:=seq(x.i=1,i=0..N);
	lprint(`dimension =`,N,`, maximal degree =`,M);
end:

# dualrep(C) returns the dual representation, and
# rank(C) returns the dimension of C

dualrep:=proc(C) sort(expand(subs(dualsubstring,C))) end:
rank   :=proc(C) sort(expand(subs(ranksubstring,C))) end:

# A monomial ideal is represented formally by the sum of its 
# monomial generators in a given degree. Hence x0+x1 is the
# (linear part of the) ideal generated by x0 and x1.

# simpleweights(A) returns the sum of one-dimensional 
# representations, one for each occurring weight in A. For
# example, simpleweights(V[j]*V[k]) should return V[j+k].
# This is useful for computing twists of a monomial ideal:
# if I sub V[d], then simpleweights(I*V[1]) is the ideal I in 
# degree d+1.

simpleweights:=proc(p) local mon;
	expand(p); if p=0 then RETURN(0) fi;
	coeffs(",variables,'mon');
	convert([mon],`+`);
end:

# twistup(J,k) returns the ideal J in degree k higher.

twistup:=proc(J,k) simpleweights(J*V[k]) end:

# Usual operations in the representation ring are taken
# care of by standard algebraic operations like +, -, *
# in conjunction with dualrep. 

Hom:=proc(A,B) sort(expand(dualrep(A)*B)) end:

# permut(elem,lis) is the action of a permutation "elem" 
# on a list "lis". The permutation is represented as a 
# product of disjoint cycles, as a list of lists. (As in
# the Maple "group" package).

permut:=proc(elem,ls)  local k,res,m,cyc,i,top;
res:=ls;
for m from 1 to nops(elem) do 
	cyc:=elem[m];
   top:=res[cyc[1]];
   for i from 1 to nops(cyc)-1 do
    res:=subsop(cyc[i]=res[cyc[i+1]],res)
   od;
  res:=subsop(cyc[i]=top,res):
od;
res
end:

# Let C^* -> T be the one-parameter subgroup given by weights 
# w0,w1,w2,w3. Thus the character xi restricts to t^(wi), i=0..3
# We need to find the products of all weights at each fixpoint.

prodwts:=proc(f) local t,cof,mon,res,i;
  cof:=[coeffs(f,variables,mon)];
  mon:=subs(x0=t^w0,x1=t^w1,x2=t^w2,x3=t^w3,[mon]);
  res:=1;
  for i from 1 to nops(mon) do 
    res:=res*subs(t=1,diff(mon[i],t))^cof[i] 
  od;
  res
end:

###############################################################

# Start by studying the fixpoints in given P^3. We go up to degree
# 9 only because we want to compute the number of twisted cubics on
# a general degree 9 hypersurface in P^7. For the Calabi-Yau three-
# folds it would suffice with M=5.

M:=9: 
setuptorus(3,M);	

# The eight types of fixpoints:

# Cohen-Macaulay types 1--4, with tangent spaces and global sections:

EE[1]:=x0*x1+x1*x2+x2*x3: 
EE[2]:=x0*x1+x1*x2+x0*x2: 
EE[3]:=x0*x1+x2^2+x0*x2:  
EE[4]:=x0^2+x0*x1+x1^2:  
for i to 4 do 
  F[i]:=expand(EE[i]*V[1]-simpleweights(EE[i]*V[1])):
  TH[i]:=expand(Hom(F[i],EE[i])*V[1]-Hom(F[i],F[i])-Hom(EE[i],EE[i])+1):
  OC[1,i]:=V[1]:
  for j from 2 to M do OC[j,i]:=V[j]-twistup(EE[i],j-2) od:
od:


# Non Cohen-Macaulay types 5--8, with tangent spaces and global sections:

beta:=expand(x0*(x0+x1+x2)):               # the same net "beta" of quadrics
q   :=[x1*x2*x3,x1*x2^2,x2^2*x3,x2^3]:     # four different cubic forms
TI  :=expand((x1+x2+x3)/x0+x3/x1+x3/x2):   # tangent space of I at beta
VF :=subs(x0=0,twistup(x1^2+x1*x2+x2^2,1)):# singular cubics through x1=x2=0

for i from 5 to 8 do
  Id[i]  :=twistup(beta,1)+q[i-4]:         # ideal in degree 3
  TH[i]  :=TI+expand(VF/q[i-4]-1+q[i-4]/(x0*x1*x2)): # (formula (x.x))
  OC[1,i]:=V[1]:
  OC[2,i]:=V[2]-beta:		              
  for j from 3 to M do OC[j,i]:=V[j]-twistup(Id[i],j-3) od:
od:         

# Isotropy groups: the permutations of the variables x0,x1,x2,x3
# that leave each T-fixpoint fixed. Permutations are products of 
# disjoint cycles, all are sugroups of the symmetric group S4.

with(group):                                # use "group" package
isot[1]:=permgroup(4,{[[1,4],[2,3]]}):      # switch x0<->x3 and x1<->x2
isot[2]:=permgroup(4,{[[1,2]],[[1,2,3]]}):  # all perms of x0,x1,x2
isot[3]:=permgroup(4,{[]}):                 # no isotropy
isot[4]:=permgroup(4,{[[1,2]],[[3,4]]}):    # switch x0<->x1 or x2<->x3
isot[5]:=permgroup(4,{[[2,3]]}):            # symmetric in x1 and x2
isot[6]:=permgroup(4,{[]}):                 # no isotropy
isot[7]:=permgroup(4,{[]}):                 # no isotropy
isot[8]:=permgroup(4,{[]}):                 # no isotropy

# To generate all fixpoints, we move the 8 given ones with one 
# representatives for each coset of the isotropy group.

S4:=permgroup(4,{[[1,2]],[[1,2,3,4]]}):     # symmetric group on x0,x1,x2,x3
for i to 8 do orb[i]:=cosets(S4,isot[i]) od:# repr. modulo isotropy 

# To prepare input in Bott's formula, here is the part of the denominator
# which comes from deformations within the fixed P^3. 

for i from 1 to 8 do den[i]:=prodwts(TH[i]) od:

# For the numerators, we need the weight products of the spaces of 
# global sections of various twists. Here they are up to degree M=9:

for i from 1 to 8 do for j from 1 to M do
      EEE[j,i]:=prodwts(OC[j,i])
od: od:

####################################################################
# Here follows the actual coding of Bott's formula. If "degrs" is a 
# list of degrees, the number of twisted cubics on the corresponding 
# complete intersection in P^N is returned. See examples at the end.

numberoftwistedcubics:=proc(N,degrs) 
local genf,WG,i,k,kk,j,fire,W4,P3,C,PN,rest,lis,result,
          TG,sstring,z,t,W,num,contP3;

# Construct "random" weight vector 
# (WG := 4, 11, 17, 32, 55, 95, 160, 267)

WG:=4,11; for i to N-1 do WG:=WG,WG[i]+WG[i+1]+2*i od;

# Generate all invariant P^3-s, i.e. all 4-element subsets
# of the variables x0,...,xN. Do this via a generating function:

genf   :=expand(product('1+z.i*t','i'=1..N+1)):
fours  :=convert(coeff(genf,t,4),list):   # products of 4 vars.
PN     :=coeff(genf,t,N+1):               # product of all vars.
sstring:=seq(z.i=i,i=1..N+1):             # vars. -> indices

# Numerator and fraction:
for i from 1 to 8 do
  num[i]:=product(EEE[degrs[k],i],k=1..nops(degrs));
  fra[i]:=unapply(num[i]/den[i],w0,w1,w2,w3);
od:

# Summation over all fixpoints in a given P3 (generic weights):

contP3:=0;
for i to 8 do for j to nops(orb[i]) do
  contP3:=contP3+fra[i](op(permut(orb[i][j],[w0,w1,w2,w3])))
od od;
contP3:=unapply(contP3,w0,w1,w2,w3);

# Main loop, over "fours", the set of 4-element subsets:

result:=0:                                # initialize result
for P3 in fours do                        # P3 represents 4 variables
   C:=PN/P3:                              # C is the remaining vars.
   fire:=subs(sstring,convert(P3,list)):  # just the indices of m
   rest:=subs(sstring,convert(C,list)):   # the remaining indices
   TG:=product(product(WG[fire[k]]-WG[rest[kk]],
                k=1..4),kk=1..N-3):       # tangent space of Grass(3,N)
   W:=seq(WG[fire[k]],k=1..4): k:='k':
   result:=result+contP3(W)/TG:           # 1/TG is a common factor
od:
result;
end:
    
numberoftwistedcubics(4,[5]);
numberoftwistedcubics(5,[3,3]);
numberoftwistedcubics(5,[4,2]);
numberoftwistedcubics(6,[3,2,2]);
numberoftwistedcubics(7,[2,2,2,2]);
numberoftwistedcubics(7,[9]);

quit;