############################################################################### # # GAP program to compute broccoli invariants with the Caporaso-Harris formula # Andreas Gathmann # December 4, 2010 # ############################################################################### # # Broccoli invariants can be computed with # Broccoli (d,r,s,alpha,beta); # where the lists alpha and beta have to be given in square brackets. For # example, # Broccoli (3,4,1,[0,1],[1]); # gives the result -8. The program is not really optimized for speed; it is # probably not very well suited high degrees of the curves. # # To compute Welschinger invariants, one can use # Welschinger (d,r,s); # which is just a short-hand notation for Broccoli (d,r,s,[],[d]). # # There are three debugging options that can be enabled / disabled # independently here: # # If you set this to 1, all the individual terms in the Caporaso-Harris formula # will be shown when computing an invariant. output := 0; # If you set this to 1, invariants with at least one real and complex point # are computed using both formulas; an error will be raised if the results # differ (which should not happen). If you set this to 0, the real equation is # used if there is at least one real point, and the complex equation otherwise. both := 0; # If you set this to 1, the algorithm is checked against some known invariants # when this file is read. checks := 0; # ############################################################################### # Functions to keep a list of the computed invariants # # To store the invariant (d,r,s,alpha,beta) as having the value n, use # SetInvariant (d,r,s,alpha,beta,n); # To get a previously computed invariant (d,r,s,alpha,beta), use # GetInvariant (d,r,s,alpha,beta); # this returns false if the number hasn't been computed yet, and the # (rational) number otherwise. # There are no restrictions on the values of d,r,s,alpha,beta except that d,r,s # must be natural numbers and alpha,beta must be lists of natural numbers. N := rec (); # Transform d,r,s,alpha,beta to a string so that it can be used to name an # invariant ToString := function (d,r,s,alpha,beta) local str; str := JoinStringsWithSeparator ([String(d),String(r),String(s),String(alpha),String(beta)],"__"); RemoveCharacters (str," []"); return ReplacedString (str,",","_"); end; GetInvariant := function (d,r,s,alpha,beta) local str; str := ToString (d,r,s,alpha,beta); if IsBound (N.(ToString (d,r,s,alpha,beta))) then return N.(ToString (d,r,s,alpha,beta)); else return false; fi; end; SetInvariant := function (d,r,s,alpha,beta,n) N.(ToString (d,r,s,alpha,beta)) := n; end; ############################################################################### # Delete trailing zeroes from a list NormalizeSequence := function (alpha) while Length(alpha)>0 and alpha[Length(alpha)]=0 do Remove (alpha); od; end; # Normalize the sequences alpha,beta and check if d,r,s,alpha,beta are valid # values so that a broccoli invariant is assigned to it. Modifies alpha,beta # and returns true/false CheckInvariant := function (d,r,s,alpha,beta) if d<=0 or r<0 or s<0 then return false; fi; if Length(alpha)>d or Length(beta)>d then return false; fi; if [1..d]*(alpha+beta+[0])<>d then return false; fi; NormalizeSequence (alpha); NormalizeSequence (beta); return r+2*s = 2*d+Sum(beta)-1; # Dimension condition end; # Return a unit vector. Input: a number i>=1, output: the list [0,0,...,0,1] # with the 1 at the i-th position UnitVector := function (i) local l; l := ShallowCopy (Zero ([1..i-1])); Add (l,1); return l; end; # Computes all possible splittings of a given degree d into alpha and beta # sequences so that I(alpha+beta)<=d. Returns a list of objects [alpha,beta] # where alpha,beta are the corresponding lists AlphaBetaList := function (d) local a,b,i,j,x,AB; AB := []; for j in [0..d] do for x in PartitionTuples (j,2) do a := []; b := []; for i in x[1] do a := a + UnitVector(i); od; for i in x[2] do b := b + UnitVector(i); od; Add (AB,[a,b]); od; od; return AB; end; # Compute a product of factorials; ProductFactorial ([k_1,...,k_r]) returns # k_1! ... k_r! ProductFactorial := function (k) return Product ([1..Length(k)], i->Factorial (k[i])); end; # Compute a multinomial coefficient; Multinomial ([k_1,...,k_r]) returns # (k_1+...+k_r)! / (k_1! ... k_r!) Multinomial := function (k) return Factorial (Sum (k)) / Product ([1..Length(k)], i->Factorial (k[i])); end; ############################################################################### # Initializes the computation of all splittings occurring in the C,...,F terms # of the Caporaso-Harris formula # Returns a record e that contains all the necessary data and can be passed to # NextSplitting # Takes the values d,r,s,alpha,beta as in the Caporaso-Harris formula. In # addition, it takes a number term in {3,4,5,6} corresponding to the terms # C,D,E,F; the splittings are then computed as needed in this term. InitSplitting := function (d,r,s,alpha,beta,term) local e; e := rec ( l := 0, dmax := d, rmax := r, smax := s, amax := alpha, bmax := beta, term := term, DIndex := 0, SIndex := 0, SList := [], kk := Length (beta)); if term=3 then e.DList := Partitions (d,2); else e.DList := Partitions (d); fi; return e; end; # Computes the next splitting (i.e. the next term in the sum) and its # coefficient # Call as NextSpltting (e); returns true/false if there is a/no further # splitting # Useable fields of the record e are: # e.l: the number of external components # e.d, e.r, e.s, e.alpha, e.beta: lists of size e.l containing the # corresponding values of the components # e.k: list of size l of the numbers k_i in the Caporaso-Harris formula # e.kk (only for term=5): k in the Caporaso-Harris formula # e.kk (only for term=6): We take the special end to be the k-th one instead of # always the first one as this is easier to implement. This makes the initial # 1/(l-1)! factor a 1/l! factor, and introduces an additional summation over k. NextSplitting := function (e) local i,c; # In this loop we try to find the next valid splitting. First advance to # the next splitting without checking if this new splitting is valid. while true do # Take the next splitting of the complex points. e.SList keeps a list # of the possible splittings; e.SIndex denotes the current position in the # list. The list is recreated below if the number e.l of external # components changes. if e.SIndex0 and e.k[i]>=e.d[i] do e.k[i] := 1; i := i-1; od; if i>0 then e.k[i] := e.k[i]+1; elif e.term=5 and e.kk0 and e.ABIndex[i]>=Length (e.ABList[i]) do e.ABIndex[i] := 1; i := i-1; od; if i>0 then e.ABIndex[i] := e.ABIndex[i]+1; else # No next splitting of the alpha and beta points available. Advance # to the next degree splitting (if there is any) and recreate the # lists of e.k, e.SList, and e.ABList. if e.DIndex1); e.ABIndex := List ([1..e.l],k->1); e.SList := OrderedPartitions (e.smax+e.l,e.l)-1; e.ABList := List ([1..e.l],k->AlphaBetaList (e.d[k])); else return false; fi; fi; # The alpha and beta splittings have been changed. Get e.alpha and # e.beta for easier access. e.alpha := List ([1..e.l], k->e.ABList[k][e.ABIndex[k]][1]); e.beta := List ([1..e.l], k->e.ABList[k][e.ABIndex[k]][2]); # The sequence alpha - sum_i alpha_i. e.asum := e.amax-Sum (e.alpha); fi; # Set e.bred to e.beta, except that we replace beta_i by beta_i - e_{k_i} # where this is required in the condition on the beta_i. e.bred := ShallowCopy (e.beta); if e.term<>3 then for i in [1..e.l] do if e.term<>6 or i<>e.kk then e.bred[i] := e.bred[i]-UnitVector(e.k[i]); fi; od; fi; # Compute the right hand side of the beta condition. e.bsum := e.bmax; if e.term=3 then e.bsum := e.bsum-UnitVector (e.k[1]+e.k[2]); fi; if e.term=5 then e.bsum := e.bsum-UnitVector (e.kk); fi; if ForAny ([1..Length (e.bsum)], k->e.bsum[k]<0) then continue; fi; # The difference of the two sides in the beta equation, must be 0 below. e.bsum := e.bsum-Sum (e.bred); fi; # The splitting of the complex points has changed. Compute the splitting of # the real points so that the dimension condition is satisfied. e.s := e.SList[e.SIndex]; e.r := List ([1..e.l], k->2*e.d[k]+Sum (e.beta[k])-1-2*e.s[k]); # A possible next splitting is computed, now check the conditions. # The degree condition for each component. This is normally # I(alpha_i+beta_i)=d_i for all i, but alpha needs to be replaced by # alpha+e_{k_i} in some cases according to the formula. c := 0; for i in [1..e.l] do if e.term=3 or (e.term=6 and i=e.kk) then if [1..e.dmax]*(e.alpha[i]+e.beta[i]+[0])+e.k[i]<>e.d[i] then c:=1; fi; else if [1..e.dmax]*(e.alpha[i]+e.beta[i]+[0])<>e.d[i] then c:=1; fi; fi; od; if c=1 then continue; fi; # The alpha condition: in term (C) alpha - sum_i alpha_i must be zero. In # all other terms it need only be non-negative. if e.term=3 then if ForAny ([1..Length(e.asum)], k->e.asum[k]<>0) then continue; fi; else if ForAny ([1..Length(e.asum)], k->e.asum[k]<0) then continue; fi; fi; # The beta condition: first check if any of the e.bred entries is negative. if ForAny ([1..e.l], k->e.k[k]<=Length (e.bred[k]) and e.bred[k][e.k[k]]<0) then continue; fi; # Now check equality of the two sides. if ForAny ([1..Length (e.bsum)], k->e.bsum[k]<>0) then continue; fi; # Check that the numbers of real and complex points on the components is OK if ForAny ([1..e.l], k->e.r[k]<0) then continue; fi; if e.rmax<>Sum (e.r) then continue; fi; # In (C), at least one of e.k[1] and e.k[2] needs to be odd if e.term=3 and not IsInt ((e.k[1]+1)*(e.k[2]+1)/2) then continue; fi; # We have found a new valid splitting. break; od; # We have a new valid splitting. Compute the factor needed in the # Caporaso-Harris formula. # First the multinomial coefficients for s and r. e.factor := Multinomial (e.s)*Multinomial (e.r); # Now the multinomial coefficients for alpha. Note that in the notation of # the paper e.max contains the sequence alpha, and e.asum contains alpha'. e.factor := e.factor * ProductFactorial (e.amax)/ProductFactorial (e.asum); e.factor := e.factor / Product ([1..e.l], i->ProductFactorial (e.alpha[i])); # The symmetry factor. We have the 1/l! which accounts for relabeling of the # components. But we only sum over _ordered_ partitions of d, so all that's # left of this 1/l! factor is the product over all 1/c_i!, where c_i denotes # the number of external components of degree i. e.factor := e.factor / Product ([1..e.dmax], i->Factorial (Number (e.d, j->j=i))); # Product over (-m)^{\alpha'_m} if e.term<>3 then e.factor := e.factor * Product ([2,4..2*Int(Length(e.asum)/2)], i->(-i)^e.asum[i]); fi; # Product over k_j and beta_{k_j}^j if e.term<>3 then for i in [1..e.l] do if e.term=6 and i=e.kk then continue; fi; if IsInt (e.k[i]/2) then e.factor := e.factor * e.k[i]; fi; e.factor := e.factor * e.beta[i][e.k[i]]; od; fi; # M_k if e.term=5 then if IsInt (e.kk/2) then e.factor := -e.factor; else e.factor := e.kk*e.factor; fi; fi; # \tilde M_{k_1} if e.term=6 then if not IsInt (e.k[e.kk]/2) then e.factor := e.k[e.kk]*e.factor; fi; fi; return true; end; ############################################################################### # This is the main recursive function to compute the broccoli invariants. # Setting debug to 2 will recompute the invariant in any case (even if has # been computed before) and output the terms used for the computation; setting # it to 1 will just print the invariant, setting it to a smaller value disables # debugging output. Compute := function (d,r,s,alpha,beta,debug) local e,k,k1,k2,n,nr,nc,m; # Print the invariant if debugging if debug=2 then Print ("\nComputing the invariant "); fi; if debug>=1 then Print ("(",d,",",r,",",s,","); n := String (alpha); RemoveCharacters (n," "); Print (n,","); n := String (beta); RemoveCharacters (n," "); Print (n,") "); fi; if debug=2 then Print ("...\n\n"); fi; # Check if the invariant has been computed previously already. if not CheckInvariant (d,r,s,alpha,beta) then Error ("Invalid invariant"); fi; if debug<2 then n := GetInvariant (d,r,s,alpha,beta); if IsRat (n) then return n; fi; fi; # Move real point, if there is any if r>0 then if debug=2 then Print ("Move real point:\n"); fi; nr := 0; # Term (A) for k in [1,3..2*Int((d-1)/2)+1] do if k<=Length(beta) and beta[k]>0 then if debug=2 then Print ("A: 1 "); fi; m := Compute (d,r-1,s,alpha+UnitVector(k),beta-UnitVector(k),debug-1); if debug=2 then Print ("-> ",m,"\n"); fi; nr := nr + m; fi; od; # Term (D) e := InitSplitting (d-1,r-1,s,alpha,beta,4); while NextSplitting (e) do if debug=2 then Print ("D: ",e.factor," "); fi; m := e.factor * Product ([1..e.l], i->Compute (e.d[i],e.r[i],e.s[i],e.alpha[i],e.beta[i],debug-1)); if debug=2 then Print ("-> ",m,"\n"); fi; nr := nr + m; od; n := nr; if debug=2 then Print ("Result from moving real point: ",nr,"\n\n"); fi; fi; # Move complex point, if there is any if (both=0 and r=0) or (both=1 and s>0) then if debug=2 then Print ("Move complex point:\n"); fi; nc := 0; # Term (B) for k in [1,3..2*Int((d-1)/2)+1] do if k<=Length(beta) and beta[k]>1 then if debug=2 then Print ("B: ",-1/2," "); fi; m := - 1/2 * Compute (d,r,s-1,alpha+UnitVector(2*k),beta-2*UnitVector(k),debug-1); if debug=2 then Print ("-> ",m,"\n"); fi; nc := nc + m; fi; od; for k1 in [1..d] do for k2 in [k1+1..d] do if k2<=Length(beta) and beta[k1]>0 and beta[k2]>0 and IsInt ((k1+1)*(k2+1)/2) then if debug=2 then Print ("B: ",-1," "); fi; m := - Compute (d,r,s-1,alpha+UnitVector(k1+k2), beta-UnitVector(k1)-UnitVector(k2),debug-1); if debug=2 then Print ("-> ",m,"\n"); fi; nc := nc + m; fi; od; od; # Term (C) e := InitSplitting (d,r,s-1,alpha,beta,3); while NextSplitting (e) do if debug=2 then Print ("C: ",e.factor," "); fi; m := e.factor * Product ([1..e.l], i->Compute (e.d[i],e.r[i],e.s[i], e.alpha[i]+UnitVector(e.k[i]),e.beta[i],debug-1)); if debug=2 then Print ("-> ",m,"\n"); fi; nc := nc + m; od; # Term (E) e := InitSplitting (d-1,r,s-1,alpha,beta,5); while NextSplitting (e) do if debug=2 then Print ("E: ",e.factor," "); fi; m := e.factor * Product ([1..e.l], i->Compute (e.d[i],e.r[i],e.s[i],e.alpha[i],e.beta[i],debug-1)); if debug=2 then Print ("-> ",m,"\n"); fi; nc := nc + m; od; # Term (F) e := InitSplitting (d-1,r,s-1,alpha,beta,6); while NextSplitting (e) do if debug=2 then Print ("F: ",e.factor," "); fi; m := e.factor; for k in [1..e.l] do if k=e.kk then m := m*Compute (e.d[k],e.r[k],e.s[k], e.alpha[k]+UnitVector(e.k[k]),e.beta[k],debug-1); else m := m*Compute (e.d[k],e.r[k],e.s[k],e.alpha[k],e.beta[k],debug-1); fi; od; if debug=2 then Print ("-> ",m,"\n"); fi; nc := nc + m; od; n := nc; if debug=2 then Print ("Result from moving complex point: ",nc,"\n\n"); fi; fi; # If the invariant has been computed both ways, check that the results agree if both=1 and r>0 and s>0 and nr<>nc then Error ("Real/complex results differ"); fi; SetInvariant (d,r,s,alpha,beta,n); return n; end; Broccoli := function (d,r,s,alpha,beta) return Compute (d,r,s,alpha,beta,2*output); end; Welschinger := function (d,r,s) return Broccoli (d,r,s,[],[d]); end; ############################################################################### # Performing initial checks with known values. if checks=0 then quit; fi; Check := function (d,r,s,alpha,beta,n) if Broccoli (d,r,s,alpha,beta) <> n then Print ("\nValue of "); Compute (d,r,s,alpha,beta,1); Print ("should be ",n,"\n"); Error ("Check failed"); fi; end; Print ("Checking some known values"); # Check all broccoli invariants for degree up to 3. Check (1,1,0,[1],[],1); Check (1,2,0,[],[1],1); Check (1,0,1,[],[1],1); Print ("."); Check (2,3,0,[2],[],1); Check (2,1,1,[2],[],1); Check (2,3,0,[0,1],[],-2); Check (2,1,1,[0,1],[],-2); Check (2,4,0,[1],[1],1); Check (2,2,1,[1],[1],1); Check (2,0,2,[1],[1],1); Check (2,5,0,[],[2],1); Check (2,3,1,[],[2],1); Check (2,1,2,[],[2],1); Check (2,4,0,[],[0,1],0); Check (2,2,1,[],[0,1],0); Check (2,0,2,[],[0,1],-1); Print ("."); Check (3,5,0,[3],[],6); Check (3,3,1,[3],[],4); Check (3,1,2,[3],[],2); Check (3,5,0,[1,1],[],-8); Check (3,3,1,[1,1],[],-4); Check (3,1,2,[1,1],[],0); Check (3,5,0,[0,0,1],[],3); Check (3,3,1,[0,0,1],[],1); Check (3,1,2,[0,0,1],[],-1); Check (3,6,0,[2],[1],8); Check (3,4,1,[2],[1],6); Check (3,2,2,[2],[1],4); Check (3,0,3,[2],[1],2); Check (3,6,0,[0,1],[1],-12); Check (3,4,1,[0,1],[1],-8); Check (3,2,2,[0,1],[1],-4); Check (3,0,3,[0,1],[1],0); Check (3,7,0,[1],[2],8); Check (3,5,1,[1],[2],6); Check (3,3,2,[1],[2],4); Check (3,1,3,[1],[2],2); Check (3,6,0,[1],[0,1],0); Check (3,4,1,[1],[0,1],0); Check (3,2,2,[1],[0,1],0); Check (3,0,3,[1],[0,1],0); Check (3,8,0,[],[3],8); Check (3,6,1,[],[3],6); Check (3,4,2,[],[3],4); Check (3,2,3,[],[3],2); Check (3,0,4,[],[3],0); Check (3,7,0,[],[1,1],0); Check (3,5,1,[],[1,1],0); Check (3,3,2,[],[1,1],0); Check (3,1,3,[],[1,1],0); Check (3,6,0,[],[0,0,1],3); Check (3,4,1,[],[0,0,1],1); Check (3,2,2,[],[0,0,1],-1); Check (3,0,3,[],[0,0,1],-3); Print ("."); # The Welschinger invariants for degree 4 Check (4,11,0,[],[4],240); Check (4,9,1,[],[4],144); Print ("."); Check (4,7,2,[],[4],80); Check (4,5,3,[],[4],40); Print ("."); Check (4,3,4,[],[4],16); Check (4,1,5,[],[4],0); Print ("."); # The Welschinger invariants for degree 5 Check (5,14,0,[],[5],18264); Print ("."); Check (5,12,1,[],[5],9096); Print ("."); Check (5,10,2,[],[5],4272); Print ("."); Check (5,8,3,[],[5],1872); Print ("."); Check (5,6,4,[],[5],744); Print ("."); Check (5,4,5,[],[5],248); Print ("."); Check (5,2,6,[],[5],64); Print ("."); Check (5,0,7,[],[5],64); Print ("."); # Two Welschinger invariants for degree 6 Check (6,17,0,[],[6],2845440); Print ("."); Check (6,15,1,[],[6],1209600); Print ("."); # The d<=3 Gromov-Witten invariants as broccoli invariants Check (2,0,2,[],[0,1],-1*2^0*1); Check (4,0,5,[],[2,1],+2*2^3*1); Print ("."); Check (6,0,8,[],[4,1],-3*2^6*12); Print ("."); # Hopefully everything was OK... Print (" OK\n");