A program in Sage which calculates numerators of Hilbert series

A link to a published Sage worksheet

Homogenous = SFAHomogeneous(QQ)
Elementary = SFAElementary(QQ)

def e(n,deg):
    return Polynomial(Elementary([deg]).expand(n)) if deg > 0 else (1 if deg == 0 else 0)
def h(n, deg):
    return Polynomial(Homogenous([deg]).expand(n)) if deg > 0 else (1 if deg == 0 else 0)
def H(n,s,l):
    return sum(map(lambda r: h(n,s-r)*e(n,r)*(-1)^r,
                             range(0,l+1)))

def a(n,k,l):
    def a_help(n,k,l,b):
        return sum(map(lambda a: ((-1)^a)*e(n,a)*H(n,k+b-a,b),
                                 range(0, k+l+1)))
    return sum(map(lambda b: (z[n]^b)*a_help(n,k,l,b), range(0,n)))

# main_hilbert calculates coefficients aij
def main_hilbert(size, n):
    if n == 2 or n==1:
        return Matrix(Polynomial, size, 1, lambda i,j: 1 if i==0 else 0)
    M = Matrix(Polynomial, size, size, lambda i,j: a(n-1, i-j, j))
    return M*main_hilbert(size, n-1)

# find_hilbert_series finds the numerators of the Hilbert series
def find_hilbert_series(n):
    v = [Matrix(Polynomial, n+1, 1, lambda i,j: z[n]^i) for n in range(0,n+1)]
    return (v[n].transpose()*main_hilbert(n+1, n))[0][0]

# Example
num_var = 5
Polynomial = PolynomialRing(QQ, num_var, 'z')
z = PolynomialRing(QQ, num_var, 'z').gens()
res = find_hilbert_series(num_var-1)



A program in Sage which calculates the Euler characteristic of invertible sheaves on spaces with four blown-up points

A link to a published Sage worksheet

num_var = 5; precision = 25
num_pairs = Integer(num_var*(num_var-1)/2)
z = PolynomialRing(QQ, num_var, 'z').gens()

# function base_change takes coefficients of a divisor \sum a_{i,j} E_{i,j}
# and returns the coefficients in the basis E_{0,1},E_{0,2},E_{0,3}, E_{0,4},E_{1,2}.
def base_change(a01,a02,a03,a04,a12,a13,a14,a23,a24,a34):
    return [a01 + a23 + a24 + a34, a02 + a13 + a14 + a34,
            a03-a13-a23-a34,
            a04-a14-a24-a34,
            a12+a13+a14+a23+a24+a34]

# kapranov_iso function takes coefficients of the divisor in basis
# E_{0,1},E_{0,2},E_{0,3}, E_{0,4},E_{1,2} and returns the gradation
# of the corresponding monomial of Grassmannian
def kapranov_iso(a01,a02,a03,a04,a12):
    return [a01+a02+a03+a04, a01, a02, a03+a12, a04+a12]

# calculation of the denominator of the Hilbert Series
denom = reduce(lambda x,y: x*y, [1 - z[i]*z[j]
                                 for i in range(num_var)
                                 for j in range(i+1,num_var)])

# fixing precision
P = PowerSeriesRing(QQ, 5, "z", default_prec = precision)
# calculation of Hilbert series
H = (P)(-z[0]^2*z[1]^2*z[2]^2*z[3]^2*z[4]^2 + z[0]^2*z[1]*z[2]*z[3]*z[4] +
         z[0]*z[1]^2*z[2]*z[3]*z[4] + z[0]*z[1]*z[2]^2*z[3]*z[4] +
         z[0]*z[1]*z[2]*z[3]^2*z[4] + z[0]*z[1]*z[2]*z[3]*z[4]^2 -
         z[0]*z[1]*z[2]*z[3] - z[0]*z[1]*z[2]*z[4] - z[0]*z[1]*z[3]*z[4] -
         z[0]*z[2]*z[3]*z[4] - z[1]*z[2]*z[3]*z[4] + 1)/denom

# calculation of the dictionary of coefficients of H
Hdict = H.dict()

# We are looking for a quadratic form in variables a01,a02,a03,a04,a12
# monomial_values returns values of subsequent monomials
# of degree less or equal to two in variables a01,a02,a03,a04,a12
def monomial_values(a):
    b = [1] + a
    return [b[i]*b[j] for i in range(num_var+1)
                      for j in range(i,num_var+1)]

# base_vector returns a vector of size num_pairs, which has a one
# at the i-th position and zeroes at all other positions
def base_vector(i):
    v = vector(QQ,num_pairs)
    v[i] = 1
    return v
# one_vector returns a vector of size num_pairs consisting only of ones
def one_vector():
    return vector(QQ, [1 for i in range(num_pairs)])

# list_of_divs returns ununique coefficients of chosen
# divisors (-2K +- E_{i,j} +- E_{k,l})
def list_of_divs():
    return [one_vector() + s1*base_vector(k) + s2*base_vector(l)
            for k in range(num_pairs)
            for l in range(k,num_pairs)
            for s1 in [1,-1]
            for s2 in [1,-1]]

# div_into_eq takes the divisor's coefficients and returns a vector
# with values of the monomials
def div_into_eq(a):
    return monomial_values(base_change(*a))

# eqmatrix is a matrix of the system of equations
eqmatrix = Matrix(map(lambda a : div_into_eq(a), list_of_divs()))

# solution is a vector of values of the quadratic form
solution = vector(map(lambda a: Hdict[
    sage.rings.polynomial.polydict.ETuple(
            kapranov_iso(*base_change(*a)))], list_of_divs()))

eqmatrix.solve_right(solution)