################################ gw #################################
# 
# This procedure computes degree d (possibly disconnected) Gromov-Witten
# invariants of a genus g curve. The syntax is as follows: gw(g,d,insert)
# computes the invariants with n insertions given by the vector
# insert=[v_1,...,v_n] where v_i is the ith insertion. v_i is either a
# number or a two vector and gives a descendant or hodge type insertion as
# follows: 
#
# To get the descendant t_k(1), set v_i=k
# To get t_k(pt) set v_i= -k (unless k=0 then use v_i=-1/2)
# To get t_k(a_i), set v_i=[k,i]
# To get t_k(b_i), set v_i=[k,-i]  
# To get ch_k(Rpi_*(f^*(O))), set v_i=[k]
# To get ch_k(Rpi_*(f^*(K))), set v_i=[-k]  (see below to change degree)
#
# (a_1,...,a_g,b_1,...,b_g) is the standard basis for H^1 of the the curve.
#
# For example to get <t_1(p)^5> we simply do gw(g,d,[-1,-1,-1,-1,-1])
# whereas gw(g,d,[4,4,[2,1],[3,-1],[-2]]) computes
# <t_4(1)t_4(1)t_2(a1)t_3(b1)ch_2(Rpi_*f^*(K))>.
#
# In order to guarantee that the sign of gw is correct, the a-loop entries
# must all preceed the b-loop entries. (Recall that the descendants of the
# loops anti-commute with each other).
#
# Note that the program computes disconnected Gromov-Witten invariants,
# i.e. integrals over the moduli space of stable maps where the domain is
# possibly disconnected. The disconnected theory is more natural from the
# point of view of the algorithm and the usual (connected) invariants are
# determined by the disconnected invariants (and vice versa). Note also
# that the euler characteristic of the domain curve is not an input to the
# procedure since it is determined by g, d, and n.
#
# The procedure gw runs as a front end for the programs GWomega,
# GWhodge, and GW. It simply reformats the inputs and calls the other
# programs.
#
# OPTIONAL DEGREE CHANGE FOR OMEGA:
#
# Optional chern insertion for a degree l bundle: If one of the
# insertions in the vector insert is [[deg]], then the insertions of
# the form [-k] will correspond to ch_k(Rpi_*(f^*(L))) where L is a
# line bundle of degree deg. The default (if there is no [[deg]]
# insertion) is deg=2g-2, i.e. L=K.

gw:=proc(g,d,insert)
local i, omega, hodge, tau_1, tau_a, tau_b, tau_p, evilsign, omega_degree;
global fast;

fast:=true;

omega_degree:=2*g-2;
evilsign:=1;
omega:=[];
hodge:=[];
tau_1:=[];
tau_a:=[];
tau_b:=[];
tau_p:=[];
for i from 1 to nops(insert) do;
  if type(insert[i],list)=false then
    if insert[i]>=0 then  
    tau_1:=[op(tau_1),insert[i]];
    else
    tau_p:=[op(tau_p),floor(-insert[i])];
    fi;
  elif nops(insert[i])=1 then
    if type(insert[i][1],list)=true then
    omega_degree:=insert[i][1][1];
    elif insert[i][1]>0 then
    hodge:=[op(hodge),insert[i][1]];
    else 
    omega:=[op(omega),-insert[i][1]];
    fi;
  elif insert[i][2]>0 then
  tau_a:=[op(tau_a),[insert[i][1],insert[i][2]]];
  else
  tau_b:=[op(tau_b),[insert[i][1],-insert[i][2]]];
  fi;
od;

if nops(tau_a)<>0 then
#WARNING("for now, gw may give the wrong sign if the a-loop entries do not all preceed the b-loop entries");
fi;

if nops(omega)<>0 then
#WARNING("running GWomega(%1,%2,%3,%4,%5,%6,%7,%8)",g,d,omega,hodge,tau_1,tau_a,tau_b,tau_p);
return GWomega(g,d,omega,hodge,tau_1,tau_a,tau_b,tau_p,omega_degree);
elif nops(hodge)<>0 then 
#WARNING("running GWhodge(%1,%2,%3,%4,%5,%6,%7)",g,d,hodge,tau_1,tau_a,tau_b,tau_p);
return GWhodge(g,d,hodge,tau_1,tau_a,tau_b,tau_p);
else
#WARNING("running GW(%1,%2,%3,%4,%5,%6)",g,d,tau_1,tau_a,tau_b,tau_p);
return GW(g,d,tau_1,tau_a,tau_b,tau_p);
fi;
end;



########################## GWomega  ##################


# GWomega computes degree d GW invariants of a genus g curve with
# arbitrary decendant insertions and also insertions of hodge classes
# and 'omega classes'.  By Hodge and omega clases we mean insertions
# of the chern character of the hodge bundle and Rpi_*(f^*(omega))
# respectively. The syntax is as follows
# GWomega(g,d,tau_omega,tau_hodge, tau_1, tau_a, tau_b, tau_p,deg) has
# the same syntax as the procedures GWhodge and GW but the additional
# vector tau_omega is a vector of numbers corresponding to the degree
# of the chern character of the omega bundle that is inserted. The
# final argument, deg, is a new addition: it is the degree of the
# bundle omega. For the original application, deg=2*g-2, i.e. omega is
# the canonical bundle, but this allows for any degree bundle.

GWomega:=proc(g,d,tau_omegain,tau_hodge, tau_1, tau_a, tau_b, tau_p,deg)
local tau_omega, tau_new, no, nh, n1, na, nb, np,
      i,  aset, bset, tlist, answer, k, ia, l,
      tau_1replace, tau_areplace, tau_breplace, tau_preplace,
      tau_1remove, tau_rem;
option remember;

#reorder tau_omega for maximum efficiency

tau_omega:=sort(tau_omegain);

no:=nops(tau_omega);
nh:=nops(tau_hodge);
n1:=nops(tau_1);
na:=nops(tau_a);
nb:=nops(tau_b);
np:=nops(tau_p);



# invariant will be 0 if a loops and b loops don't pair up. We also
# check that the inputs are consistant.

#if na<>nb then 
#WARNING("Number of a loops and b loops are different"); return 0 
#fi;
#  for i from 1 to na do;
#    if nops(tau_a[i])<> 2 or nops(tau_b[i])<>2
#    then ERROR(" %1 and %2 must each have 2 entries",tau_a[i],tau_b[i]);
#    fi;
#    if type(g,integer) then
#      if (tau_a[i][2]>g or tau_b[i][2]>g) then
#      ERROR("Loop index out of range. Neither %1 nor %2 should not be greater than g=%3", tau_a[i][2],tau_b[i][2],g); 
#      fi;
#    fi;
#  od;
#aset:={seq(tau_a[i][2],i=1..na)};
#bset:={seq(tau_b[i][2],i=1..nb)};
#  if not(aset subset bset and bset subset aset) then
#  WARNING("should be automatically 0");
#  return 0;
#  fi;

# set invariant to 0 if a tau subscript is negative. Negative
# subscripts will occur in the recursion

tlist:={op(tau_1),seq(tau_a[i][1],i=1..na),seq(tau_b[i][1],i=1..nb),op(tau_p)};
   if member(-1,tlist) 
   then return 0 
   fi;

# if there are no Omega insertions, then use GWhodge to evaluate

if no=0 then return GWhodge(g,d,tau_hodge, tau_1,tau_a,tau_b,tau_p) fi;


# Apply the operator H_{k} to remove the first omega entry. There are
# three cases H_1, H_{2k-1}, H_{2k}.

tau_new:=[op(2..no,tau_omega)];

if tau_omega[1]=1 then
# apply H_1
# WARNING("applying H_1}");

answer:=-deg/2*
        GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_b,[1,op(tau_p)],deg);

  for i from 1 to n1 do;
  tau_1remove:=rementry(tau_1,i);
  l:=tau_1[i];
  answer:=answer + deg/2*
      GWomega(g,d,tau_new,tau_hodge,tau_1remove,tau_a,tau_b,[l,op(tau_p)],deg);
  od;

answer:=answer + 1/12*
        GWomega(g,d,tau_new,tau_hodge,[2,op(tau_1)],tau_a,tau_b,tau_p,deg);

  for i from 1 to n1 do;
  tau_1replace:=tau_1;
  tau_1replace[i]:=tau_1[i]+1;
  answer:=answer - 1/12*
          GWomega(g,d,tau_new,tau_hodge,tau_1replace,tau_a,tau_b,tau_p,deg);
  od; 

  for i from 1 to np do;
  tau_preplace:=tau_p;
  tau_preplace[i]:=tau_p[i]+1;
  answer:=answer - 1/12*
          GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_b,tau_preplace,deg);
  od; 

  for i from 1 to na do;
  tau_areplace:=tau_a;
  tau_breplace:=tau_b;
  tau_areplace[i][1]:=tau_a[i][1]+1;
  tau_breplace[i][1]:=tau_b[i][1]+1;
  answer:=answer - 1/12*
          GWomega(g,d,tau_new,tau_hodge,tau_1,tau_areplace,tau_b,tau_p,deg);
  answer:=answer - 1/12*
          GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_breplace,tau_p,deg);
  od;

answer:=answer+ 1/12*
    GWomega(g,d,tau_new,tau_hodge,[0,op(tau_1)],tau_a,tau_b,[0,op(tau_p)],deg);

  for i from 1 to g do;
  answer:=answer + 1/12*(-1)^(na+1)*
          GWomega(g,d,tau_new,tau_hodge,tau_1,
                  [[0,i],op(tau_a)],[[0,i],op(tau_b)],tau_p,deg);
  od;

  for i from 1 to no-1 do;
    if tau_new[i]=1 then
    tau_rem:=rementry(tau_new,i);
    answer:=answer- 1/24*deg*
        GWomega(g,d,tau_rem,tau_hodge,tau_1,tau_a,tau_b,[0,0,op(tau_p)],deg);  
    fi;
  od;

fi;

if tau_omega[1]>1 and type(tau_omega[1],odd) then

# apply H_{2k-1}

k:=(tau_omega[1]+1)/2;
# WARNING("applying H_{2k-1} with k=%1",k);

answer:=GWomega(g,d,tau_new,tau_hodge,[2*k,op(tau_1)],tau_a,tau_b,tau_p,deg);

  for i from 1 to n1 do;
  tau_1replace:=tau_1;
  tau_1replace[i]:=tau_1[i]+2*k-1;
  answer:=answer-
          GWomega(g,d,tau_new,tau_hodge,tau_1replace,tau_a,tau_b,tau_p,deg);
  od;

  for i from 1 to np do;
  tau_preplace:=tau_p;
  tau_preplace[i]:=tau_p[i]+2*k-1;
  answer:=answer-
          GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_b,tau_preplace,deg);
  od;

  for i from 1 to na do;
  tau_areplace:=tau_a;
  tau_areplace[i][1]:=tau_a[i][1]+2*k-1;
  tau_breplace:=tau_b;
  tau_breplace[i][1]:=tau_b[i][1]+2*k-1;
  answer:=answer-
          GWomega(g,d,tau_new,tau_hodge,tau_1,tau_areplace,tau_b,tau_p,deg)-
          GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_breplace,tau_p,deg);
  od;

  for l from 0 to 2*k-2 do;
  answer:=answer + (-1)^l*
  GWomega(g,d,tau_new,tau_hodge,
          [l,op(tau_1)],tau_a,tau_b,[2*k-2-l,op(tau_p)],deg);
    for i from 1 to g do;
    answer:= answer - (-1)^(l+na)*
    GWomega(g,d,tau_new,tau_hodge,
            tau_1,[[l,i],op(tau_a)],[[2*k-2-l,i],op(tau_b)],tau_p,deg);
    od;
  od;

  for i from 1 to no-1 do;
#print("foo");
    if tau_new[i]=1 then
    tau_rem:=rementry(tau_new,i);
      for l from 0 to 2*k-2 do;
      answer:=answer-deg/2*(-1)^l*
      GWomega(g,d,tau_rem,tau_hodge,
                  tau_1,tau_a,tau_b,[l,2*k-2-l,op(tau_p)],deg);  
      od;
    fi;
  od;

answer:=bernoulli(2*k)/(2*k)!*answer;
fi;

if tau_omega[1]>1 and type(tau_omega[1],even) then
# apply H_{2k}

k:=tau_omega[1]/2;
# WARNING("applying H_{2k} with k=%1",k);

answer:=GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_b,[2*k,op(tau_p)],deg);

  for i from 1 to n1 do
  tau_1remove:=rementry(tau_1,i);
  l:=tau_1[i];
  answer:=answer -
   GWomega(g,d,tau_new,tau_hodge,
               tau_1remove,tau_a,tau_b,[l+2*k-1,op(tau_p)],deg);
  od;

  for l from 0 to 2*k-2 do;
  answer:=answer + 1/2*(-1)^l*
   GWomega(g,d,tau_new,tau_hodge,tau_1,tau_a,tau_b,[l,2*k-2-l,op(tau_p)],deg);
  od;


answer:=deg*bernoulli(2*k)/(2*k)!*answer;
fi;

if type(g,integer)=false then answer:=factor(answer) fi;

return answer;
end;


########################## GWhodge  ##################

# GWhodge computes degree d GW invariants of a genus g curve with
# arbitrary decendant insertions and also insertions of hodge classes,
# namely insertions of the (2i-1)th chern character of the hodge
# bundle. The syntax is as follows GWhodge(g,d,tau_hodge, tau_1,
# tau_a, tau_b, tau_p) has the same syntax as the procedure GW but the
# additional vector tau_hodge is a vector of odd numbers corresponding
# to the degree of the chern character of the Hodge bundle that is
# inserted.

GWhodge:=proc(g,d,tau_hodgein, tau_1in, tau_a, tau_b, tau_p)
local tau_new, tau_hodge, tau_1, nh, n1, na, nb, np,
      i, isodd, aset, bset, tlist, answer, k, GWhsum, ia,
      l, tau_1replace, tau_areplace, tau_breplace, tau_preplace;
option remember;

# reorder entries to maximize efficiency
tau_hodge:=sort(tau_hodgein);
tau_1:=sort(tau_1in);

nh:=nops(tau_hodge);
n1:=nops(tau_1);
na:=nops(tau_a);
nb:=nops(tau_b);
np:=nops(tau_p);

# check that the tau_hodge entries are odd

isodd:=proc(n) type(n,odd);end;
  if nops(remove(isodd,tau_hodge))<>0 then 
#  WARNING("one of the hodge insertions is not an odd number"); 
  return 0;
  fi;

# invariant will be 0 if a loops and b loops don't pair up. We also
# check that the inputs are consistant.

#if na<>nb then 
#WARNING("Number of a loops and b loops are different"); return 0 
#fi;
#  for i from 1 to na do;
#    if nops(tau_a[i])<> 2 or nops(tau_b[i])<>2
#    then ERROR("%1 and %2 must each have 2 entries",tau_a[i],tau_b[i]);
#    fi;
#    if type(g,integer) then
#      if tau_a[i][2]>g or tau_b[i][2]>g then
#      ERROR("Loop index out of range. Neither %1 nor %2 should not be greater than g=%3", tau_a[i][2],tau_b[i][2],g) 
#      fi;
#    fi;
#  od;
#aset:={seq(tau_a[i][2],i=1..na)};
#bset:={seq(tau_b[i][2],i=1..nb)};
#  if not(aset subset bset and bset subset aset) then
#  WARNING("should be automatically 0");
#  return 0;
#  fi;

# set invariant to 0 if a tau subscript is negative. Negative
# subscripts will occur in the recursion

tlist:={op(tau_1),seq(tau_a[i][1],i=1..na),seq(tau_b[i][1],i=1..nb),op(tau_p)};
   if member(-1,tlist) 
   then return 0 
   fi;

# if there are no Hodge insertions, then use GW to evaluate

if nh=0 then return GW(g,d,tau_1,tau_a,tau_b,tau_p) fi;

# Apply the operator D_{2k-1} to remove the first hodge entry

k:=(tau_hodge[1]+1)/2;
tau_new:=[op(2..nh,tau_hodge)];

# single derivative terms

answer:=GWhodge(g,d,tau_new,[2*k,op(tau_1)],tau_a,tau_b,tau_p);

for i from 1 to n1 do;
tau_1replace:=tau_1;
tau_1replace[i]:=tau_1[i]+(2*k-1);
answer:=answer-GWhodge(g,d,tau_new,tau_1replace,tau_a,tau_b,tau_p);
od;

for i from 1 to np do;
tau_preplace:=tau_p;
tau_preplace[i]:=tau_p[i]+(2*k-1);
answer:=answer-GWhodge(g,d,tau_new,tau_1,tau_a,tau_b,tau_preplace);
od;

for i from 1 to na do;
tau_areplace:=tau_a;
tau_areplace[i][1]:=tau_a[i][1]+2*k-1;
answer:=answer - GWhodge(g,d,tau_new,tau_1,tau_areplace,tau_b,tau_p);
od;

for i from 1 to nb do;
tau_breplace:=tau_b;
tau_breplace[i][1]:=tau_b[i][1]+2*k-1;
answer:=answer - GWhodge(g,d,tau_new,tau_1,tau_a,tau_breplace,tau_p);
od;

# Double derivative terms. Sign in odd term is to commute the
# tau_a-tau_b insertion from being next to each other, to being in my
# standard format.

for l from 0 to 2*k-2 do;
GWhsum:=GWhodge(g,d,tau_new,[op(tau_1),l],tau_a,tau_b,[op(tau_p),2*k-2-l]);
  for ia from 1 to g do;
  GWhsum:=GWhsum - (-1)^na*
  GWhodge(g,d,tau_new,tau_1,[[l,ia],op(tau_a)],[[2*k-2-l,ia],op(tau_b)],tau_p);
  od;
answer:=answer + (-1)^l*GWhsum;
od;

return bernoulli(2*k)/(2*k)!*answer;
end;

##########################  GW  #########################
#
# GW computes the GW invariants of a target curve of genus g with
# arbitrary insertions of decendants of points, loops, and
# constants. It uses the (conjectural) Virasoro relations to remove
# decendants of constants and express the invariant in terms of
# invariants with only loop and point decendants. This is then
# evaluated with GWloop which uses the (proven) formulas of Rahul and Andrei.

# Syntax is as follows. GW(g,d,tau_1,tau_a,tau_b,tau_p) returns the GW
# invariant of degree d, target genus g, point decendant insertions
# given by the vector tau_p, loop insersions given by the vectors
# tau_a and tau_b, and the decendants of 1 given by the vector
# tau_1. The vectors tau_p and tau_1 are vectors of non-negative
# integers and the vectors tau_a and tau_b are vectors whose entries
# are two vectors whose first entry is the decendant number and whose
# second number indicates which a-loop or b-loop is used. For example,
# GW(g,d,[1,2],[[3,4],[5,6]],[[7,8],[9,10]],[11]) evaluates
# <t_1(1)t_2(1)t_3(a4)t_5(a6)t_7(b8)t_9(b10)t_11(p)>_d with target
# genus g where a1,...,ag,b1,...,bg is a standard basis for H_1 of the
# target.

GW:=proc(g,d,tau_1in,tau_a,tau_b,tau_p) 
local tau_1, n1, na, nb, np, aset, bset, i,k, answer, tlist,
      j, tau_new, tau_1shift, tau_ashift, tau_bshift, tau_pshift,
      zeroindex1, zeroindexp, numzero1, numzerop, tau_1remove, tau_premove,
      ia, ib, tau_aremove, tau_bremove, evilsign, multiplier, zeros,
      tau_newreplace, tau_areplace, tau_breplace, tau_preplace,l;
option remember;

#reorder tau_1 for maximum efficiency
tau_1:=sort(tau_1in);

n1:=nops(tau_1);
na:=nops(tau_a);
nb:=nops(tau_b);
np:=nops(tau_p);


# invariant will be 0 if a loops and b loops don't pair up. We also
# check that the inputs are consistant.

#if na<>nb then 
#WARNING("Number of a loops and b loops are different"); return 0 
#fi;
#  for i from 1 to na do;
#    if nops(tau_a[i])<> 2 or nops(tau_b[i])<>2
#    then ERROR("%1 and %2 must each have 2 entries",tau_a[i],tau_b[i]);
#    fi;
#    if tau_a[i][2]>g or tau_b[i][2]>g then
#    ERROR("Neither %1 nor %2 should not be greater than %3",
#    tau_a[i][2],tau_b[i][2],g) 
#    fi;
#  od;
#aset:={seq(tau_a[i][2],i=1..na)};
#bset:={seq(tau_b[i][2],i=1..nb)};
#  if not(aset subset bset and bset subset aset) then
#  WARNING("should be automatically 0");
#  return 0;
#  fi;

# set invariant to 0 if a tau subscript is negative. Negative
# subscripts will occur in the recursion

tlist:={op(tau_1),seq(tau_a[i][1],i=1..na),seq(tau_b[i][1],i=1..nb),op(tau_p)};
   if member(-1,tlist) 
   then return 0 
   fi;

# if there are no tau_1 terms then use GWloop to evaluate otherwise
# use the L_k recursion operator to remove a tau_1 term

if n1=0 then 
  if fast=true then
  return fastGWloop(g,d,np,tau_p,na,tau_a,tau_b)
  else return GWloop(g,d,np,tau_p,na,tau_a,tau_b)
  fi;
fi;

k:=tau_1[1]-1;
tau_new:=[op(2..n1,tau_1)];

# now apply the appropriate L_k.  The three main cases are k=-1, k=0,
# and k>0.

if k=-1 then 
#WARNING("applying L_{-1} with tau_new=%1,tau_p=%2",tau_new,tau_p);

answer:=0;

# derivative terms of L_{-1}:

  for j from 1 to n1-1 do;
  tau_1shift:=tau_new;
  tau_1shift[j]:=tau_1shift[j]-1;
  answer:=answer+ GW(g,d,tau_1shift,tau_a,tau_b,tau_p);
  od;

  for j from 1 to na do;
  tau_ashift:=tau_a;
  tau_ashift[j][1]:=tau_ashift[j][1]-1;
  answer:=answer+GW(g,d,tau_new,tau_ashift,tau_b,tau_p);
  od;

  for j from 1 to nb do;
  tau_bshift:=tau_b;
  tau_bshift[j][1]:=tau_bshift[j][1]-1;
  answer:=answer+GW(g,d,tau_new,tau_a,tau_bshift,tau_p);
  od;

  for j from 1 to np do;
  tau_pshift:=tau_p;
  tau_pshift[j]:=tau_pshift[j]-1;
  answer:=answer+GW(g,d,tau_new,tau_a,tau_b,tau_pshift);
  od;

# constant terms of L_{-1}

zeroindex1:=0;
zeroindexp:=0;
numzero1:=0;
numzerop:=0;
  for j from 1 to n1-1 do;
    if tau_new[j]=0 then
    numzero1:=numzero1+1;
    zeroindex1:=j;
    fi;
  od;
  for j from 1 to np do;
    if tau_p[j]=0 then
    numzerop:=numzerop+1; 
    zeroindexp:=j;
    fi;
  od;

  if numzero1<>0 and numzerop<>0 then
  tau_1remove:=rementry(tau_new,zeroindex1);
  tau_premove:=rementry(tau_p,zeroindexp);
  answer:=answer+
          numzero1*numzerop*GW(g,d,tau_1remove,tau_a,tau_b,tau_premove);
  fi;

# odd constant term contributions of L_{-1}

  for ia from 1 to na do;
    if tau_a[ia][1]=0 then
    for ib from 1 to nb do;
      if tau_b[ib][1]=0 and tau_b[ib][2]=tau_a[ia][2] then
      tau_aremove:=rementry(tau_a,ia);
      tau_bremove:=rementry(tau_b,ib);
      evilsign:=(-1)^(ia-1+na+ib);
      answer:=answer+evilsign*GW(g,d,tau_new,tau_aremove,tau_bremove,tau_p);
      fi;
    od;
    fi;
  od;

fi;


if k=0 then
#WARNING("applying L_0 with tau_new=%1, tau_a=%2, tau_b=%3, tau_p=%4",tau_new,tau_a,tau_b,tau_p);

answer:= (2*g-2)*GW(g,d,tau_new,tau_a,tau_b,[0,op(tau_p)]);
               
multiplier:= sum('tau_new[ii]','ii'=1..n1-1) +
             sum('tau_a[ii][1]+1','ii'=1..na)+ 
             sum('tau_b[ii][1]','ii'=1..nb) +
             sum('tau_p[ii]+1','ii'=1..np);
             
answer:= answer + multiplier*GW(g,d,tau_new,tau_a,tau_b,tau_p);

  for i from 1 to n1-1 do;
  tau_1remove:=rementry(tau_new,i);
  answer:=answer+
          (2-2*g)*GW(g,d,tau_1remove,tau_a,tau_b,[op(tau_p),tau_new[i]-1]);
  od;

# constant term of L_0


zeros:=kinvect(tau_new,0);
  if zeros[1]>1 then
#  WARNING("doing constant term of L_0");  
  tau_1remove:=rementry(rementry(tau_new,zeros[2][2]),zeros[2][1]);
  answer:=answer+(2-2*g)*binom(zeros[1],2)*
          GW(g,d,tau_1remove,tau_a,tau_b,tau_p);
  fi;

fi;


# L_k for k>0 we divide by (k+1)! and write the products in terms of
# binomial coefficients

if k>0 then

#WARNING("Applying L_%1 with tau_new=%2, tau_a=%3, tau_b=%4, tau_p=%5",k,tau_new,tau_a,tau_b,tau_p);

answer:= (2*g-2)*sum('1/rr','rr'=1..k+1)
        *GW(g,d,tau_new,tau_a,tau_b,[k,op(tau_p)]);

  for i from 1 to n1-1 do;
  l:=tau_new[i];
  tau_newreplace:=tau_new;
  tau_newreplace[i]:=l+k;
  answer:=answer + binom(k+l,k+1)*GW(g,d,tau_newreplace,tau_a,tau_b,tau_p);
  od;

  for i from 1 to np do;
  l:=tau_p[i];
  tau_preplace:=tau_p;
  tau_preplace[i]:=l+k;
  answer:=answer+binom(k+l+1,k+1)*GW(g,d,tau_new,tau_a,tau_b,tau_preplace);
  od;

  for i from 1 to na do;
  l:=tau_a[i][1];
  tau_areplace:=tau_a;
  tau_areplace[i][1]:=l+k;
  answer:=answer + binom(k+l,k+1)*GW(g,d,tau_new,tau_areplace,tau_b,tau_p);
  od;

  for i from 1 to nb do;
  l:=tau_b[i][1];
  tau_breplace:=tau_b;
  tau_breplace[i][1]:=l+k;
  answer:=answer+binom(k+l+1,k+1)*GW(g,d,tau_new,tau_a,tau_breplace,tau_p);
  od;

  for i from 1 to n1-1 do;
  l:=tau_new[i];
  tau_1remove:=rementry(tau_new,i);
    if l=0 then
    multiplier:=1/(k+1);
    else multiplier:=binom(k+l,k+1)*sum('1/rr','rr'=l..k+l) 
    fi;
  answer:=answer+(2-2*g)*multiplier
                 *GW(g,d,tau_1remove,tau_a,tau_b,[op(tau_p),k+l-1]);
  od;

  for l from 0 to k-2 do;
  answer:=answer+ (1-g)/(k+1)/binom(k,l+1)
                   *GW(g,d,tau_new,tau_a,tau_b,[op(tau_p),l,k-l-2]);
  od;

fi;


return factor(answer);
end;

# kinvect takes in a vector vect and a number k and returns 
# [number of k's in vect,[index of 1st k,...,index of last k]]

kinvect:=proc(vect,k)
local n, j, numk, position;
n:=nops(vect);
numk:=0;
position:=[];
  for j from 1 to n do;
    if vect[j]=k then 
    numk:=numk+1;
    position:=[op(position),j];
    fi;
  od;
return [numk,position];
end;

# rementry is a procedure that removes vect[index] from index

rementry:=proc(vect,index);
if index<1 or index>nops(vect)
then ERROR("index=%1 out of range in procedure rementry",index)
fi;
[op(1..index-1,vect),op(index+1..nops(vect),vect)];
end;



#######################  GWloop  #########################
#
# reduces the GW invariants of a curve of genus g involving
# descendants of point and loop classes to GW invariants only
# involving descendants of the point.  Uses the formula of Rahul and
# Andrei.

# A basis for loop classes on the genus g surface will be given by
# a[1],..,a[g],b[1],...,b[g] with a[i]*b[j]= delta[i,j] and a[i]*a[j]=
# b[n]*b[m]=0. The only non-zero GW invariants involve an equal number
# of a and b insertions (but they could have repetition).  

# The argumets of GWloop are g the target genus, d the degree, n the
# number point insertions, m the number of loop insertions, and the
# descendant data tau_point, tau_a and tau_b. tau_point is a length n
# vector where tau_point[i] is the multiplicity of the ith point
# descendant. tau_a and tau_b are m vectors of 2 vectors so tau_a[i]
# is data for the ith a-loop insertions; it is a 2-vector [x,k]
# indicating the insertion is tau[x]a[k]. Similarly for tau_b.

# In summary, the GW that GWloop evaluates is the formula from Rahul's
# notes so for example
# GWloop(g,d,2,[2,3],3,[[1,2],[3,4],[5,6]],[[7,8],[9,10],[11,12]])
# evalutates the invariant
# <t_2(p)t_3(p)t_1(a2)t_3(a4)t_5(a6)t_7(b8)t_9(b10)t_11(b12)>_d for
# target a genus g curve


with(combinat);with(group);

GWloop:=proc(g,d,n,tau_point,m,tau_a,tau_b) 

local perms, answer, sigma, s, bin, point, cup, i, k, 
      prod_term, xi, ysigmai, new_tau, aset, bset, evilsign;
option remember;


#
# report arguments
#
#WARNING
#("running GWloop with g=%2, d=%1, n=%3, m=%4, 
#tau_point=%5, tau_a=%6, tau_b=%7",d,g,n,m,tau_point,tau_a,tau_b);
#
# check that input data is consistant
#
#  if nops(tau_point)<> n 
#  then ERROR("You haven't specified %1 point descendants",n); 
#  fi;
#  if nops(tau_a)<> m or nops(tau_b)<> m 
#  then ERROR("You haven't specified %1 loop descendants",m)
#  else 
#    for k from 1 to m do;
#      if nops(tau_a[k])<> 2 or nops(tau_b[k])<>2
#      then ERROR
#      ("%1 and %2 must each have 2 entries",tau_a[k],tau_b[k]);
#      fi;
#      if tau_a[k][2]>g or tau_b[k][2]>g then 
#      ERROR("Neither %1 nor %2 should not be greater than %3",
#      tau_a[k][2],tau_b[k][2],g) 
#      fi;
#    od; 
#  fi;
#
# if m=0 then use GWpoint to evaluate
#
  if m=0 then  return GWpoint(g,d,n,tau_point); 
  else answer:=0;
  fi;

# if the a-loop and b-loop insertions don't match up it should be
# automatically 0

#aset:={seq(tau_a[i][2],i=1..m)};
#bset:={seq(tau_b[i][2],i=1..m)};
#  if not(aset subset bset and bset subset aset) then
#  WARNING("should be automatically 0");
#  return 0;
#  fi;

# if m>8 then 
# WARNING("GWloop is being called with m=%1. Evilsign will be wrong",m);
# fi;

perms:=permute(m);

  for s from 1 to m! do;
  sigma:=perms[s];
  prod_term:=1;
    for i from 1 to m do; 
    xi:=tau_a[i][1]; 
    ysigmai:=tau_b[sigma[i]][1];
    bin:=binom(xi+ysigmai,xi);
      if tau_a[i][2]=tau_b[sigma[i]][2] then 
      cup:=1 else 
      cup:=0 
      fi;
    new_tau:=[seq(tau_a[j][1]+tau_b[sigma[j]][1]-1,j=1..m),op(tau_point)];
    prod_term:=prod_term*bin*cup;
    od;
#  evilsign:=(-1)^(m*(m-1)/2)*parity(convert(perms[s],'disjcyc'));  
# the folling formula for evilsign works for m less than 7
# evilsign:=(-1)^(m*(m-1)/2)*(-1)^(s*(s-1)/2+floor((s-1)/24));
    if prod_term<>0 then
    answer:=answer + evils(m,s)*prod_term*GWpoint(g,d,n+m,new_tau);
    fi;
  od;
answer;
end;  

binom:=proc(n,k)
option remember;
binomial(n,k);
end;

######### evils ########
#
# This procedure computes the evil sign in GWloop. The "real" formula
# for this sign is
#
# (-1)^(m*(m-1)/2)*parity(convert(permute(m)[s],'disjcyc'));
#
# but instead we use a "closed" form for this sign based on the order
# that Maple lists the set of permutations in its procedure
# permutations(m) (which is called by GWloop). Note that currently,
# this procedure only gives the right answer for m up to 8, but it is
# easy to ammend.

evils:=proc(m,s)
local sig;
option remember;

if m<5 then
sig:=(-1)^(m*(m-1)/2)*(-1)^(s*(s-1)/2);
elif m<7 then
sig:=(-1)^(m*(m-1)/2 + s*(s-1)/2 + floor((s-1)/24));
elif m<9 then
sig:=(-1)^(m*(m-1)/2 + s*(s-1)/2 + floor((s-1)/24) +floor((s-1)/720));
else ERROR("%1 is 9 or bigger",m)
fi;
sig;
end;

#######################  fastGWloop  #########################
#
# reduces the GW invariants of a curve of genus g involving
# descendants of point and loop classes to GW invariants only
# involving descendants of the point.  Uses the formula of Rahul and
# Andrei.

# A basis for loop classes on the genus g surface will be given by
# a[1],..,a[g],b[1],...,b[g] with a[i]*b[j]= delta[i,j] and a[i]*a[j]=
# b[n]*b[m]=0. The only non-zero GW invariants involve an equal number
# of a and b insertions (but they could have repetition).  

# The argumets of GWloop are g the target genus, d the degree, n the
# number point insertions, m the number of loop insertions, and the
# descendant data tau_point, tau_a and tau_b. tau_point is a length n
# vector where tau_point[i] is the multiplicity of the ith point
# descendant. tau_a and tau_b are m vectors of 2 vectors so tau_a[i]
# is data for the ith a-loop insertions; it is a 2-vector [x,k]
# indicating the insertion is tau[x]a[k]. Similarly for tau_b.

# In summary, the GW that GWloop evaluates is the formula from Rahul's
# notes so for example
# GWloop(g,d,2,[2,3],3,[[1,2],[3,4],[5,6]],[[7,8],[9,10],[11,12]])
# evalutates the invariant
# <t_2(p)t_3(p)t_1(a2)t_3(a4)t_5(a6)t_7(b8)t_9(b10)t_11(b12)>_d for
# target a genus g curve.

# This is the fast version of GWloop. It begins by reordering the
# descendant insertions so that all the insertions corresponding to
# the same loop are grouped together. It then checks to see if there
# is any redundancy of descendant number within a singling group (and
# returns zero if so). It then sums over the permutations of
# individual groupsings. This is the main place where the program
# saves time. It does not sum over the full permutation group the way
# that GWloop does.


fastGWloop:=proc(g,d,n,tau_point,m,tau_a,tau_b) 

local perms, answer, sig, S, bin, i, k, alist, blist, curag, curbg,
       x, y, new_tau, aset, bset, evilsign, reordersign;
option remember;

  if m=0 then  return GWpoint(g,d,n,tau_point); 
  else answer:=0;
  fi;

alist:=[seq([],i=1..g)];
blist:=[seq([],i=1..g)];
reordersign:=1;

for i from 1 to m do;
curag:=tau_a[i][2];
curbg:=tau_b[i][2];
reordersign:=reordersign*
            (-1)^( sum('nops(alist[ii])','ii'=curag..g)
                  +sum('nops(blist[ii])','ii'=curbg..g));
alist[curag]:=[tau_a[i][1],op(alist[curag])];
blist[curbg]:=[tau_b[i][1],op(blist[curbg])];
od;

for i from 1 to g do;
if nops(alist[i])<>nops({op(alist[i])}) then return 0 fi;
if nops(blist[i])<>nops({op(blist[i])}) then return 0 fi;
od;



perms:=[seq(permute(nops(alist[ii])),ii=1..g)];


S:=cartprod(perms);
  while not S[finished] do 
  sig:=S[nextvalue]();

  new_tau:=tau_point;
  bin:=1;
  evilsign:=1;
    for i from 1 to g do;
    evilsign:=evilsign*parity(convert(sig[i],'disjcyc'));
      for k from 1 to nops(alist[i]) do;
      x:=alist[i][k];
      y:=blist[i][sig[i][k]];
      bin:=bin*binom(x+y,x);
      new_tau:=[x+y-1,op(new_tau)];
      od;
    od;
  answer:=answer+evilsign*bin*GWpoint(g,d,nops(new_tau),new_tau);
  end do;

answer:=answer*(-1)^(m*(m-1)/2)*reordersign;
end;  




################################ GWpoint ##########################
#
# GWpoint computes the degree d GW invariants with target a genus g curve
# with n insertions which are descendants of the point class. tau[1],...,tau[n]
# is a vector of integers which gives the multiplicities of the descendants.
# Uses the formula of Rahul and Andrei

with(combinat);

GWpoint:=proc(g,d,n,tau)
#global numcalls;
local parts, nparts, Answer,j,i,prod;
option remember;

#numcalls:=numcalls+1;

#  if nops(tau)<> n 
#  then ERROR("You haven't specified %1 descendants",n); 
#  fi;


parts:=partition(d);
nparts:=nops(parts);
Answer:=0;
  for j from 1 to nparts do;
  prod:=1;
    for i from 1 to n do;
    prod:=prod*power(parts[j],tau[i]+1)/(tau[i]+1)!;
    od;
  Answer:=Answer+(d!/Chi(parts[j],parts[1]))^(2*g-2)*prod;
  od;
end;

################# power ###############

power:=proc(part,l)
local c,pow,length,i;
option remember;

# Note that the partitions are passed with there entries in increasing
# order but the formula for power is more natural with the partitions
# in decreasing order. This is handled in this ruitine by using the
# negative scripting: part[-i] gives the ith entry *from the end of part*.



length:=nops(part);
c:=coeftayl(convert(series(z/2/sinh(z/2),z=0,l+4),polynom),z=0,l+1);
pow:=0;
  for i from 1 to length do;
  pow:=pow+(part[-i]-i+1/2)^l-(-i+1/2)^l;
  od;
pow:=pow+l!*c;

end;




############## branch ################

# The procedures branch2, branch4,...,branchX compute the degree d GW
# invariant of a local genus g curve for branch number X. In other
# words it is evaluating the chern class c_X of the bundle
# Rpi_*f^*(L+L) where L is a degree g-1 line bundle. It does it by
# expressing this chern class in terms of chern character insertions
# which are then evaluated by the procedure gw. These procedures are
# called by the procedures Ztil and Z to produce the partition
# functions (see below). The actual formula in branchX is produced by
# the procedure chernnum included at the bottom.

branch2:=proc(g,d)
local AN;
AN:=      gw(g,d,[[-2],[[g-1]]])
    +     gw(g,d,[[-2],[[g-1]]])
    + 1/2*gw(g,d,[[-1],[-1],[[g-1]]])
    + 1/2*gw(g,d,[[-1],[-1],[[g-1]]])
    +     gw(g,d,[[-1],[-1],[[g-1]]]);
end;

branch4:=proc(g,d)
local AN;
AN:=   (1/24)* gw(g,d,[[-1],[-1],[-1],[-1],[[g-1]]])      
      +(1/6 )* gw(g,d,[[-1],[-1],[-1],[-1],[[g-1]]])    
      +(1/4 )* gw(g,d,[[-1],[-1],[-1],[-1],[[g-1]]])
      +(1/6 )* gw(g,d,[[-1],[-1],[-1],[-1],[[g-1]]])
      +(1/24)* gw(g,d,[[-1],[-1],[-1],[-1],[[g-1]]])
      +(1/2 )* gw(g,d,[[-1],[-1],[-2],[[g-1]]])
      +(1/2 )* gw(g,d,[[-1],[-1],[-2],[[g-1]]])
      +(1   )* gw(g,d,[[-1],[-1],[-2],[[g-1]]])
      +(1   )* gw(g,d,[[-1],[-1],[-2],[[g-1]]])
      +(1/2 )* gw(g,d,[[-1],[-1],[-2],[[g-1]]])
      +(1/2 )* gw(g,d,[[-1],[-1],[-2],[[g-1]]])
      +(2   )* gw(g,d,[[-1],[-3],[[g-1]]])
      +(2   )* gw(g,d,[[-1],[-3],[[g-1]]])
      +(2   )* gw(g,d,[[-1],[-3],[[g-1]]])
      +(2   )* gw(g,d,[[-1],[-3],[[g-1]]])
      +(1/2 )* gw(g,d,[[-2],[-2],[[g-1]]])
      +(1   )* gw(g,d,[[-2],[-2],[[g-1]]])
      +(1/2 )* gw(g,d,[[-2],[-2],[[g-1]]])
      +(6   )* gw(g,d,[[-4],[[g-1]]])
      +(6   )* gw(g,d,[[-4],[[g-1]]]);
return AN;
end;

branch6:=proc(g,d)
local AN;

AN:=
           4/45 *gw(g, d, [[-1], [-1], [-1], [-1], [-1], [-1], [[g - 1]]])
         + 4/3 *gw(g, d, [[-1], [-1], [-1], [-1], [-2], [[g - 1]]])
         + 4 *gw(g, d, [[-1], [-1], [-2], [-2], [[g - 1]]])
         + 4/3 *gw(g, d, [[-2], [-2], [-2], [[g - 1]]])
         + 16/3 *gw(g, d, [[-1], [-1], [-1], [-3], [[g - 1]]])
         + 16 *gw(g, d, [[-1], [-2], [-3], [[g - 1]]])
         + 8 *gw(g, d, [[-3], [-3], [[g - 1]]])
         + 24 *gw(g, d, [[-1], [-1], [-4], [[g - 1]]])
         + 24 *gw(g, d, [[-2], [-4], [[g - 1]]])
         + 96 *gw(g, d, [[-1], [-5], [[g - 1]]])
         + 240 *gw(g, d, [[-6], [[g - 1]]]);

return AN;
end;

branch8:=proc(g,d)
local AN;

AN:=
    2/315*gw(g,d,[[-1], [-1], [-1], [-1], [-1], [-1], [-1], [-1], [[g-1]]])
   +8/45*gw(g,d,[[-1], [-1], [-1], [-1], [-1], [-1], [-2], [[g-1]]])
   +4/3*gw(g,d,[[-1], [-1], [-1], [-1], [-2], [-2], [[g-1]]])
   +8/3*gw(g,d,[[-1], [-1], [-2], [-2], [-2], [[g-1]]])
   +2/3*gw(g,d,[[-2], [-2], [-2], [-2], [[g-1]]])
   +16/15*gw(g,d,[[-1], [-1], [-1], [-1], [-1], [-3], [[g-1]]])
   +32/3*gw(g,d,[[-1], [-1], [-1], [-2], [-3], [[g-1]]])
   +16*gw(g,d,[[-1], [-2], [-2], [-3], [[g-1]]])
   +16*gw(g,d,[[-1], [-1], [-3], [-3], [[g-1]]])
   +16*gw(g,d,[[-2], [-3], [-3], [[g-1]]])
   +8*gw(g,d,[[-1], [-1], [-1], [-1], [-4], [[g-1]]])
   +48*gw(g,d,[[-1], [-1], [-2], [-4], [[g-1]]])
   +24*gw(g,d,[[-2], [-2], [-4], [[g-1]]])
   +96*gw(g,d,[[-1], [-3], [-4], [[g-1]]])
   +72*gw(g,d,[[-4], [-4], [[g-1]]])
   +64*gw(g,d,[[-1], [-1], [-1], [-5], [[g-1]]])
   +192*gw(g,d,[[-1], [-2], [-5], [[g-1]]])
   +192*gw(g,d,[[-3], [-5], [[g-1]]])
   +480*gw(g,d,[[-1], [-1], [-6], [[g-1]]])
   +480*gw(g,d,[[-2], [-6], [[g-1]]])
   +2880*gw(g,d,[[-1], [-7], [[g-1]]])
   +10080*gw(g,d,[[-8], [[g-1]]]);

return AN;
end;


branch10:=proc(g,d)
local AN;

AN:=
    4/15*gw(g,d,[[-2], [-2], [-2], [-2], [-2], [[g-1]]])
   +4/14175*gw(g,d,[[-1],[-1],[-1],[-1],[-1],[-1],[-1],[-1],[-1],[-1],[[g-1]]])
   +4/315*gw(g,d,[[-1],[-1], [-1],[-1],[-1],[-1],[-1],[-1],[-2],[[g-1]]])
   +8/45*gw(g,d,[[-1], [-1], [-1], [-1], [-1], [-1], [-2], [-2], [[g-1]]])
   +8/9*gw(g,d,[[-1], [-1], [-1], [-1], [-2], [-2], [-2], [[g-1]]])
   +4/3*gw(g,d,[[-1], [-1], [-2], [-2], [-2], [-2], [[g-1]]])
   +32/315*gw(g,d,[[-1],[-1],[-1],[-1], [-1], [-1], [-1], [-3], [[g-1]]])
   +32/15*gw(g,d,[[-1], [-1], [-1],[-1], [-1], [-2], [-3], [[g-1]]])
   +32/3*gw(g,d,[[-1], [-1], [-1], [-2], [-2], [-3], [[g-1]]])
   +32/3*gw(g,d,[[-1], [-2], [-2], [-2], [-3], [[g-1]]])
   +16/3*gw(g,d,[[-1], [-1], [-1], [-1], [-3], [-3], [[g-1]]])
   +32*gw(g,d,[[-1],[-1], [-2], [-3], [-3], [[g-1]]])
   +16*gw(g,d,[[-2], [-2], [-3], [-3], [[g-1]]])
   +64/3*gw(g,d,[[-1], [-3], [-3], [-3], [[g-1]]])
   +16/15*gw(g,d,[[-1], [-1], [-1], [-1], [-1], [-1], [-4], [[g-1]]])
   +16*gw(g,d,[[-1], [-1], [-1], [-1], [-2], [-4], [[g-1]]])
   +48*gw(g,d,[[-1], [-1], [-2], [-2], [-4], [[g-1]]])
   +16*gw(g,d,[[-2], [-2], [-2], [-4], [[g-1]]])
   +64*gw(g,d,[[-1], [-1], [-1], [-3], [-4], [[g-1]]])
   +192*gw(g,d,[[-1], [-2], [-3], [-4], [[g-1]]])
   +96*gw(g,d,[[-3], [-3], [-4], [[g-1]]])
   +144*gw(g,d,[[-1], [-1], [-4], [-4],[[g-1]]])
   +144*gw(g,d,[[-2], [-4], [-4], [[g-1]]])
   +64/5*gw(g,d,[[-1], [-1], [-1], [-1], [-1], [-5], [[g-1]]])
   +128*gw(g,d,[[-1], [-1], [-1], [-2], [-5], [[g-1]]])
   +192*gw(g,d,[[-1], [-2], [-2], [-5], [[g-1]]])
   +384*gw(g,d,[[-1], [-1], [-3], [-5], [[g-1]]])
   +384*gw(g,d,[[-2], [-3], [-5], [[g-1]]])
   +1152*gw(g,d,[[-1], [-4], [-5], [[g-1]]])
   +1152*gw(g,d,[[-5], [-5], [[g-1]]])
   +160*gw(g,d,[[-1], [-1], [-1], [-1], [-6], [[g-1]]])
   +960*gw(g,d,[[-1], [-1], [-2], [-6], [[g-1]]])
   +480*gw(g,d,[[-2], [-2], [-6], [[g-1]]])
   +1920*gw(g,d,[[-1], [-3], [-6], [[g-1]]])
   +2880*gw(g,d,[[-4], [-6], [[g-1]]])
   +1920*gw(g,d,[[-1], [-1], [-1], [-7], [[g-1]]])
   +5760*gw(g,d,[[-1], [-2],[-7], [[g-1]]])
   +5760*gw(g,d,[[-3], [-7], [[g-1]]])
   +20160*gw(g,d,[[-1], [-1], [-8], [[g-1]]])
   +20160*gw(g,d,[[-2], [-8], [[g-1]]])
   +161280*gw(g,d,[[-1], [-9], [[g-1]]])
   +725760*gw(g,d,[[-10], [[g-1]]]);

return AN;
end;




######## N_0  ########
#
# computes N_0^b(g) using Rahul and Carel's formula for the integral
# of lambda_{h-1}^3 over M_h which appears in the virtual class
# calculation of the degree 0 local invariant Z_0^b(g). Note that the
# "branch number" b is defined by b=(2h-2)-d(2g-2) where h is the
# genus of the domain curve.

N_0:=proc(g,b);
if type(b,even)=false then ERROR("b=%1 must be even",b) fi;
return (-1)^(b/2+1)*(3*g-3)*abs(bernoulli(b+2)*bernoulli(b))/(b+2)/b/b!;
end;


####### Ztil #####

# computes the series Z~_d(g) to order N. The ~ indicates that it
# includes degree zero components.

Ztil:=proc(g,d,N)
local AN;

AN:= gw(g,d,[]);

if N>0 then
AN:=AN+branch2(g,d)*t^2;
fi;

if N>2 then
AN:=AN+branch4(g,d)*t^4;
fi;

if N>4 then
AN:=AN+branch6(g,d)*t^6;
fi;

if N>6 then
AN:=AN+branch8(g,d)*t^8;
fi;

if N>8 then
AN:=AN+branch10(g,d)*t^10;
fi;

if N>10 then
ERROR("No can do more than branch number 10");
fi;

return AN;
end;

######## Z #######

# Z computes the partition function Z_d(g) of the degree d theory of a
# local genus g curve in a CY3 to order t^N. This theory does *not*
# have maps with degree 0 components. It is obtained from the full
# theory Z~ by removing the constant contributions.

Z:=proc(g,d,N)
local N0;
N0:=sum('-N_0(g,2*k)*t^(2*k)','k'=1..N+1);
series(exp(N0)*Ztil(g,d,N),t=0,N+2);
end;



################# chernnum ########################

# this procedure produces the formula for evaluating the nth chern
# class of Rpi_*f^*(L+L) in terms of the chern character insertions
# for the degree g-1 line bundle L. It is used to produce the code in
# the procedures branch2, branch4, etc...

chernnum:=proc(n)
local parts, j, monosub, P, i, k, vect, mono;

parts:=partition(n);

monosub:={};
for j from 1 to nops(parts) do;
P:=[seq(0,l=1..n)];
  for k from 1 to nops(parts[j]) do;
  P[parts[j][k]]:=P[parts[j][k]]+1;
  od;

vect:=[];
mono:=1;
  for i from 1 to n do;
  mono:=mono*O[i]^(P[i]);
  vect:=[op(vect),seq([-i],l=1..P[i])];
  od;
monosub:={op(monosub),mono=GWsub(g,d,[op(vect),[[g-1]]])};

od;


subs(monosub,
simplify(expand(
sum('
  1/permfact(parts[p])*
    product('
    (-1)^(parts[p][k]-1)*(parts[p][k]-1)!*(-O[parts[p][k]]-O[parts[p][k]])
    ',
    'k'=1..nops(parts[p]))  
  ',
'p'=1..nops(parts))
)));

end;


permfact:=proc(a)
local d,k,i,A;
option remember;

d:=sum(a[i],i=1..nops(a));
   for k from 1 to d do;
   A[k]:=0;
      for i from 1 to nops(a) do;
      if a[i]=k then A[k]:=A[k]+1 end if;
      od;
   od;
product(A[j]!,j=1..d);
end;

############# Y  ##############

# the procedure Y takes the Z series, puts in the renomalizing factor
# of t^(d(2g-2)) and then expresses it in terms of the variable
# y=-4sin^2(t/2).

Y:=proc(g,d,N);
series(subs(t=2*arcsin(1/2*sqrt(-y)),t^(2*d*(g-1))*Z(g,d,N)),y,d*(g-1)+N+1);
end;