Chern numbers of K3^[n] in Sage

chern_numbers.sage

NTHREADS = 8 # use this to specify the number of threads

class TnVariety:
    """
    A variety with a torus action for computing with Bott's formula
    """
    def __init__(self, n, points, weight):
        self.n = n
        self.points = points
        self.weight = weight
    @parallel(NTHREADS)
    def compute(self, points):
        ans = {p: 0 for p in self.pp}
        for pt in points:
            w = self.weight(pt)
            cherns = [1] + [chern(k, w) for k in range(1, len(w)+1)]
            ctop = cherns[-1]
            for p in ans.keys():
                ans[p] += prod(cherns[k] for k in p) * QQ((1, ctop))
        return ans
    def chern_numbers(self):
        self.pp = Partitions(self.n)
        # no parallel
        # return self.compute(tuple(self.points))
        each = max(len(self.points) // NTHREADS, 10)
        num = len(self.points) // each
        c = self.compute([self.points[i*each:(i+1)*each] for i in range(num)] + [self.points[num*each:]])
        ans = {p: 0 for p in self.pp}
        for (_, ci) in c:
            for p in ans.keys():
                ans[p] += ci[p]
        return ans

    def __mul__(self, other):
        n = self.n + other.n
        points = cartesian_product([self.points, other.points])
        weight = lambda p: self.weight(p[0]) + other.weight(p[1])
        return TnVariety(n, points, weight)

def chern(k, w):
    def dfs(k, n):
        if k < 1:
            ans[0] += pp[0]
        else:
            for m in range(k, n+1):
                pp[k-1] = pp[k] * w[m-1]
                dfs(k-1, m-1)
    ans, pp = [0], [0]*k + [1]
    dfs(k, len(w))
    return ans[0]

def Pn(n):
    """
    Projective space of dimension n, with a torus action
    """
    return TnVariety(n, range(n+1), lambda k: [x-k for x in range(0, k)] + [x-k for x in range(k+1, n+1)])

def _hilb(n, parts, wt):
    points = [x for pp in parts for x in cartesian_product([Partitions(x) for x in pp])]
    def weight(pt):
        w = []
        for k,l,m in wt:
            if len(pt[k]) > 0:
                b = pt[k].conjugate()
                for s in range(len(pt[k])):
                    j = pt[k][s]
                    for i in range(j):
                        w.append(l*(i-j)   + m*(b[i]-s-1))
                        w.append(l*(j-i-1) + m*(s-b[i]))
        return w
    return TnVariety(2*n, points, weight)

def hilb_P2(n):
    """
    Hilbert scheme of n points on the projective plane, with a torus action
    """
    return _hilb(n, [[a, b-a-1, n+1-b] for a,b in Combinations(n+2, 2)], [[0,1,n+1], [1,n,-1], [2,-n-1,-n]])

def hilb_P1xP1(n):
    """
    Hilbert scheme of n points on the product of two projective lines, with a
    torus action
    """
    return _hilb(n, [[a, b-a-1, c-b-1, n+2-c] for a,b,c in Combinations(n+3, 3)], [[0,-n,-1], [1,n,-1], [2,-n,1], [3,n,1]])

@cached_function
def _product(m, n):
    R = PolynomialRing(QQ, ["c"+str(i) for i in range(m)] + ["d"+str(i) for i in range(n)], order=TermOrder("wdeglex", list(range(1,m+1)) + list(range(1,n+1))))
    TX = 1+sum(R.gens()[0:m])
    TY = 1+sum(R.gens()[m:m+n])
    cTXY = TX * TY
    c = by_degree(cTXY, m+n)
    pp = Partitions(m+n)
    ans = {p: {} for p in pp}
    for p in pp:
        cp = prod(c[i] for i in p)
        monoms, coeffs = cp.monomials(), cp.coefficients()
        for i in range(len(monoms)):
            d = monoms[i].exponents()[0]
            if sum(d[i] * (i+1) for i in range(m)) == m:
                p1 = Partition([i+1 for i in range(m-1,-1,-1) for j in range(d[i])])
                p2 = Partition([i+1 for i in range(n-1,-1,-1) for j in range(d[i+m])])
                ans[p][(p1,p2)] = coeffs[i]
    return ans

def chern_nb_product(X, Y):
    m, n = X[0], Y[0]
    if m > n:
        return chern_nb_product(Y, X)
    d = _product(m, n)
    pp = Partitions(m+n)
    ans = {p: 0 for p in pp}
    for p in pp:
        for p1p2 in d[p].keys():
            p1, p2 = p1p2
            a = X[1][p1] if p1 in X[1].keys() else 0
            b = Y[1][p2] if p2 in Y[1].keys() else 0
            ans[p] += d[p][p1p2] * a * b
    return m+n, ans

def by_degree(x, n):
    c = [0] * (n+1)
    if x.parent().ngens() == 1:
        g = x.parent().gen()
        l = x.list()
        for i in range(max(n, len(l))):
            c[i] = l[i] * g^i
    monoms, coeffs = x.monomials(), x.coefficients()
    for i in range(len(monoms)):
        d = monoms[i].degree() 
        if d <= n:
            if type(coeffs) == dict:
                c[d] += coeffs[monoms[i]] * monoms[i]
            else:
                c[d] += coeffs[i] * monoms[i]
    return c

def capped_log(x, n):
    e = by_degree(x, n)
    p = [e[1]] + [0] * (n-1)
    for i in range(n-1):
        p[i+1] = (i+2) * e[i+2] - sum(e[j+1] * p[i-j] for j in range(i+1))
    return sum(QQ((1, i+1))*p[i] for i in range(n))

def capped_exp(x, n):
    comps = by_degree(x, n)
    p = [i * comps[i] for i in range(n+1)]
    e = [1] + [0] * (n)
    for i in range(n):
        e[i+1] = QQ((1, i+1)) * sum(p[j+1] * e[i-j] for j in range(i+1))
    return sum(e)

def hilb_K3(n):
    """
    Compute the non-trivial Chern numbers of the Hilbert scheme of n points on
    a K3 surface
    """
    S = PowerSeriesRing(QQ, ["a"+str(i+1) for i in range(n)]+["b"+str(i+1) for i in range(n)], order=TermOrder("wdeglex", tuple([i for i in range(1,n+1)]+[i for i in range(1,n+1)])), default_prec=n+1)
    A, B = gens(S)[:n], gens(S)[n:]
    HS = capped_exp(-16*capped_log(1+sum(A), n)+18*capped_log(1+sum(B), n), n)
    monoms, coeffs = HS.monomials(), HS.coefficients()
    pp = [p for p in Partitions(2*n) if all(map(is_even, p))]
    ans = {p: 0 for p in pp}
    P2n    = [hilb_P2(k+1).chern_numbers()    for k in range(n)]
    P1xP1n = [hilb_P1xP1(k+1).chern_numbers() for k in range(n)]
    for i in range(len(monoms)):
        if monoms[i].degree() == n:
            d = monoms[i].exponents()[0]
            lP2n    = [(2*k+2, P2n[k])    for k in range(n) for i in range(d[k])]
            lP1xP1n = [(2*k+2, P1xP1n[k]) for k in range(n) for i in range(d[k+n])]
            c = reduce(chern_nb_product, lP2n + lP1xP1n)[1]
            for p in pp:
                ans[p] += coeffs[monoms[i]] * c[p]
    return ans

# examples
print((Pn(1)*Pn(3)).chern_numbers()) # Chern numbers of a product
print(hilb_P2(2).chern_numbers())    # Chern numbers of P2^[2]
print(hilb_K3(3))                    # Chern numbers of K3^[3]

SageCell