GRRcomp.py

# -*- coding: utf-8 -*-
r"""
Chern character of the derived pushforward of a line bundle `\O(D)` on
the universal curve C_{g,n} over the space Mbar_{g,n} of stable curves via
the Grothendieck-Riemann-Roch (GRR) formula.

The algorithm implemented is the formula of Theorem 1 from the paper [PRvZ20]_ of
Pagani, Ricolfi and van Zelm. The formula itself is in the function
``generalized_hodge_chern(l,d,a,dmax,g,n)`` where using the notation from [PRvZ20]_:

- ``l`` is the integer l from equation (0.1)
- ``d`` is a list [d_1,...,d_n]
- ``a`` is a list of triples [h,S,a_{h,S}], where h is an integer and S a list
  so that [h,S] indicates the graph with two vertices and one edge with genus h
  on the left and markings S on the left. The integer a_{h,S} is the
  corresponding value to this graph from equation (0.1). Only the graphs with
  nonzero a_{h,S} need to be entered.

- ``dmax`` is the maximal codim to which we want to compute the formula of Theorem 1.
- ``g`` and ``n`` are respectively the genus and the number of marked points

As described in the paper [PRvZ20]_ this can be used to compute the DR cycle in
certain cases (namely when (d=[0,...,\pm 1,...,\mp 1,...,0]). The DR cycle can
be computed directly in this case by

    DR_phi(g,d)

The DR cycle can also be computed using [JPPZ17]_. This method has already been
implemented in admcycles.sage as DR_cycle(g,d). To verify that the methods give
the same result modulo Pixton relations one can type:

    (DR_phi(g,d) - DR_cycle(g,d)).is_zero()

We verified that these classes are the same for d=[1,-1] and g=1,2,3,4. For
higher values of g the computation takes too long. Note that in the case of
g=4 these are codimension 4 classes. The ring R^4(M_4,2) is generated by 3990
distinct decorated stratum classes, so this is a serious check on the
correctness of these formula's. Note also that the two methods do not give the
same result if we do not mod out by Pixton relations. i.e.

    (DR_phi(g,d) - DR_cycle(g,d)).simplify()

is not zero.
"""
import operator

from .admcycles import tautclass, decstratum, stgraph, psicl, kappacl, psiclass, kappaclass
from .admcycles import chern_char_to_poly, chern_char_to_class, fundclass

from sage.rings.all import ZZ, QQ
from sage.matrix.constructor import matrix
from sage.misc.misc import powerset
from sage.misc.misc_c import prod
from sage.arith.all import bernoulli, binomial, multinomial
from sage.combinat.all import bernoulli_polynomial, IntegerVectors


def divided_bernoulli(n):
    """
    Return Bernoulli number  of index ``n`` divided by ``n!``.

    EXAMPLES::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: divided_bernoulli(12)
        -691/1307674368000
    """
    return bernoulli(n) / ZZ(n).factorial()


def divided_bernoulli_polynomial(l, n):
    """
    Return Bernoulli polynomial of index ``n`` at ``l``, divided by ``n!``.

    EXAMPLES::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: x = polygen(QQ,'x')
        sage: divided_bernoulli_polynomial(x,4)
        1/24*x^4 - 1/12*x^3 + 1/24*x^2 - 1/720
    """
    return bernoulli_polynomial(l, n) / ZZ(n).factorial()


def handle_g_n(g, n):
    """
    Helper function that fetches the global ``g`` and ``n`` parameters.

    See also :func:`admcycles.reset_g_n`.

    The global ``g`` and ``n`` should rather be replaced by the creation
    of an object ``MMbar(g, n)`` with attributes ``g`` and ``n``.
    """
    if g is None:
        from .admcycles import g
    if n is None:
        from .admcycles import n
    return g, n


def generalized_hodge_chern(l, d, a, dmax=None, g=None, n=None):
    r"""
    Computes the Chern character of the derived pushforward of a line bundle \O(D) on
    the universal curve C_{g,n} over the space Mbar_{g,n} of stable curves.

    A divisor D on C_{g,n} up to pullbacks from Mbar_{g,n} takes the form

      D = l \tilde{K} + sum_{i=1}^n d_i \sigma_i  +  \sum_{h,S} a_{h,S} C_{h,S}

    where the numbers l, d_i and a_{h,S} are integers, \tilde{K} is the relative canonical
    class of the morphism C_{g,n} -> Mbar_{g,n}, \sigma_i is the image of the ith section
    and C_{h,S} is the boundary divisor of C_{g,n} where the moving point lies on a genus h
    component with markings given by the set S.
    For such a divisor this function computes the Chern character

      ch(R^\bullet \pi_* \O(D))

    up to codimension dmax using the formula given in [PRvZ20]_.

    NOTE: This formula assumes that the number n of marked points is at least 1.

    INPUT:

    l : integer
      the coefficient in front of \tilde{K}
    d : list
      vector of integers (d_1,...,d_n)
    a : list
      list of triples [h,S,a_{h,S}] where a_{h,S} is nonzero (and S is a subset of [1,...,n] )
    dmax : integer
      maximum codimension, set to 3g+3-n by default
    g : integer
      genus
    n : integer
      number of markings, for the formula to be correct we need at least one

    EXAMPLES::

      sage: from admcycles import *
      sage: from admcycles.GRRcomp import *
      sage: g=1;n=2
      sage: l=0;d=[1,-1];a=[[1,[],-1]]
      sage: generalized_hodge_chern(l,d,a,1,g,n)
      Graph :      [1] [[1, 2]] []
      Polynomial : 1/12*(kappa_1)_0 - 13/12*psi_1 - 1/12*psi_2
      <BLANKLINE>
      Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
      Polynomial : 1/24
      <BLANKLINE>
      Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
      Polynomial : -11/12
    """
    g, n = handle_g_n(g, n)
    if dmax is None:
        dmax = 3 * g - 3 + n

    return sum(generalized_hodge_chern_single(t, l, d, a, g, n) for t in range(1, dmax + 1))


def generalized_hodge_chern_single(t, l, d, a, g=None, n=None):
    r"""
    Computes the degree t part of the Chern character of the derived pushforward of a line
    bundle \O(D) on the universal curve C_{g,n} over the space Mbar_{g,n} of stable curves.

    A divisor D on C_{g,n} up to pullbacks from Mbar_{g,n} takes the form

      D = l \tilde{K} + sum_{i=1}^n d_i \sigma_i  +  \sum_{h,S} a_{h,S} C_{h,S}

    where the numbers l, d_i and a_{h,S} are integers, \tilde{K} is the relative canonical
    class of the morphism C_{g,n} -> Mbar_{g,n}, \sigma_i is the image of the ith section
    and C_{h,S} is the boundary divisor of C_{g,n} where the moving point lies on a genus h
    component with markings given by the set S.
    For such a divisor this function computes the degree t part

      ch_t(R^\bullet \pi_* \O(D))

    of the Chern character of R^\bullet \pi_* \O(D) using the formula given in
    [PRvZ20]_.

    NOTE: This formula assumes that the number n of marked points is at least 1.

    EXAMPLES::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: l=0; d=[1,-1]; a=[[1,[],-1]]
        sage: generalized_hodge_chern_single(1, l, d, a, 1, 2)
        Graph :      [1] [[1, 2]] []
        Polynomial : 1/12*(kappa_1)_0 - 13/12*psi_1 - 1/12*psi_2
        <BLANKLINE>
        Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
        Polynomial : 1/24
        <BLANKLINE>
        Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
        Polynomial : -11/12

    TESTS::

        sage: reset_g_n(1, 2)
        doctest:...: DeprecationWarning: reset_g_n is deprecated. Please use TautologicalRing instead.
        See https://gitlab.com/modulispaces/admcycles/-/merge_requests/109 for details.
        sage: generalized_hodge_chern_single(1, l, d, a)
        Graph :      [1] [[1, 2]] []
        Polynomial : 1/12*(kappa_1)_0 - 13/12*psi_1 - 1/12*psi_2
        <BLANKLINE>
        Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
        Polynomial : 1/24
        <BLANKLINE>
        Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
        Polynomial : -11/12
    """
    g, n = handle_g_n(g, n)
    if n == 0:
        raise ValueError('There needs to be at least one marked point for generalized_hodge_chern_single to work')

    a.sort(key=operator.itemgetter(1, 2))  # Order a to our convention
    result = generalized_hodge_chern_K(t, l, a, g, n) + generalized_hodge_chern_S(t, l,
                                                                                  d, a, g, n) + todd_inv_exp(a, t, g, n)  # -generalized_hodge_chern_exp(t,a,g,n)
    result = result.degree_part(t)
    result.simplify()
    return result


def todd_inv_exp(a, t, g=None, n=None):
    """
    EXAMPLES::

      sage: from admcycles import *
      sage: from admcycles.GRRcomp import *
      sage: todd_inv_exp([],1,1,2)
      Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
      Polynomial : 1/24
      <BLANKLINE>
      Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
      Polynomial : 1/12
    """
    g, n = handle_g_n(g, n)
    a.sort(key=operator.itemgetter(1, 2))  # Order a to our convention
    result = tautclass([], g, n)
    for t1 in range(1, t, 2):
        coeff_t1 = divided_bernoulli(t1 + 1) / ZZ(t - t1).factorial()
        for tvect in IntegerVectors(t - t1, len(a)):
            result += coeff_t1 * prod(ai[2]**tvecti for ai, tvecti in zip(a, tvect)) * \
                multinomial(tvect) * todd_inv_exp_graphs(a, t1, tvect, g, n)
    result += todd_inv_sigma(t, g, n)
    return result


##
# The case without the graphs
def todd_inv_sigma(t, g, n):
    if t % 2 == 0:
        return tautclass([], g, n)

    if n > 0:
        seprange = [[h, [1] + s] for h in range(g + 1) for s in powerset(range(2, n + 1))
                    if not (h == 0 and not s) and not (h == g and len(s) > n - 3)]
    else:
        seprange = [[h, []] for h in range(1, g // 2 + 1)]
    # The list [h,S] of possible nonisomorphic irreducible codimension 1 graphs
    decsum = sum((-psicl(n + 1, n + 2))**i * psicl(n + 2, n + 2)**(t - 1 - i) for i in range(t))  # The decoration
    tautclass1 = tautclass([decstratum(stgraph([g - 1], [list(range(1, n + 3))],
                           [(n + 1, n + 2)]), poly=decsum)], g, n)  # The irreducible graph
    tautclass2 = sum(tautclass([decstratum(stgraph([ahS[0], g - ahS[0]], [ahS[1] + [n + 1], list(set(range(1, n + 1)) - set(
        ahS[1])) + [n + 2]], [(n + 1, n + 2)]), poly=decsum)], g, n) for ahS in seprange)  # The separating graphs

    result = divided_bernoulli(t + 1) * (QQ((1, 2)) * tautclass1 + tautclass2)
    result = result.degree_part(t)
    result.simplify()
    return result


##
# The graphs in the sigma part
def todd_inv_exp_graphs(a, t1, tvect, g, n):
    r"""
    TESTS:

    We check that https://gitlab.com/modulispaces/admcycles/-/merge_requests/149 does not reappear::

        sage: from admcycles.GRRcomp import todd_inv_exp_graphs
        sage: len(todd_inv_exp_graphs([[1, [1], -1]], 1, [1], 3, 2)._terms) == 5
        True
    """
    # assume a is ordered
    tvn = [tvi for tvi in tvect if tvi != 0]
    atn = [ai for ai, tvi in zip(a, tvect) if tvi != 0]
    L = len(tvn)
    result = 0 * fundclass(g, n)
    if all(atn[i][0] <= atn[i + 1][0] and atn[i][1] <= atn[i + 1][1]
           for i in range(L - 1)):
        # restrict to the case nonzero elements of tvect give a ordered sequence in a (i.e. each element smaller than or equal to the next).

        Snc = list(powerset(set(range(1, n + 1)) - set(atn[-1][1])))  # Powerset of complement of Sn

        # The sum of different graphs (bigger than or equal (hn,Sn))
        for h in range(atn[-1][0], g + 1):
            for S in Snc:
                # the strictly bigger ones:
                if not (h == atn[-1][0] and not S) and not (h == g and n - len(atn[-1][1]) - len(S) < 2):
                    vertices = [atn[0][0]] + [atn[i][0] - atn[i - 1][0]
                                              for i in range(1, L)] + [h - atn[-1][0]] + [g - h]
                    labels = [atn[0][1] + [n + 1]] + [list(set(atn[i][1]) - set(atn[i - 1][1])) + [n + 2 * i, n + 2 * i + 1] for i in range(1, L)] + [list(
                        set(S) - set(atn[-1][1])) + [n + 2 * L, n + 2 * L + 1]] + [list(set(range(1, n + 1)) - set(atn[-1][1]) - set(S)) + [n + 2 * L + 2]]
                    involution = [(n + 2 * i + 1, n + 2 * i + 2) for i in range(L + 1)]
                    graph = stgraph(vertices, labels, involution)
                    psi = prod([sum([binomial(tvn[i] - 1, j) * (-1)**(tvn[i] - 1) * psicl(n + 2 * i + 1, n + 2 * L + 2)**j * psicl(n + 2 * i + 2, n + 2 * L + 2)**(tvn[i] - 1 - j) for j in range(tvn[i])]) for i in range(L)] + [sum(
                        [(-psicl(n + 2 * L + 1, n + 2 * L + 2))**j * psicl(n + 2 * L + 2, n + 2 * L + 2)**(t1 - 1 - j) for j in range(t1)])])  # psi class of graph coinciding graphs. Note last one in a is only on one side and to the power of a[-1].
                    result += tautclass([decstratum(graph, poly=psi)], g, n)

                # the case where they are equal
                if (h == atn[-1][0] and not S):
                    vertices = [atn[0][0]] + [atn[i][0] - atn[i - 1][0] for i in range(1, L)] + [g - atn[-1][0]]
                    labels = [atn[0][1] + [n + 1]] + [list(set(atn[i][1]) - set(atn[i - 1][1])) + [n + 2 * i, n + 2 * i + 1]
                                                      for i in range(1, L)] + [list(set(range(1, n + 1)) - set(atn[-1][1])) + [n + 2 * L]]
                    involution = [(n + 2 * i + 1, n + 2 * i + 2) for i in range(L)]
                    graph = stgraph(vertices, labels, involution)
                    psi = prod([sum([binomial(tvn[i] - 1, j) * (-1)**(tvn[i] - 1) * psicl(n + 2 * i + 1, n + 2 * L)**j * psicl(n + 2 * i + 2, n + 2 * L)**(tvn[i] - 1 - j) for j in range(tvn[i])])
                               for i in range(L - 1)] + [(-psicl(n + 2 * L - 1, n + 2 * L))**(tvn[-1])] + [sum([(-psicl(n + 2 * L - 1, n + 2 * L))**j * psicl(n + 2 * L, n + 2 * L)**(t1 - 1 - j) for j in range(t1)])])
                    result += tautclass([decstratum(graph, poly=psi)], g, n)

        # The case of the irreducible graph in Sigma:
        if g - atn[-1][0] > 0:
            vertices = [atn[0][0]] + [atn[i][0] - atn[i - 1][0] for i in range(1, L)] + [g - atn[-1][0] - 1]
            labels = [atn[0][1] + [n + 1]] + [list(set(atn[i][1]) - set(atn[i - 1][1])) + [n + 2 * i, n + 2 * i + 1] for i in range(
                1, L)] + [list(set(range(1, n + 1)) - set(atn[-1][1])) + [n + 2 * L, n + 2 * L + 1, n + 2 * L + 2]]
            involution = [(n + 2 * i + 1, n + 2 * i + 2) for i in range(L + 1)]
            graph = stgraph(vertices, labels, involution)
            psi = prod([sum([binomial(tvn[i] - 1, j) * (-1)**(tvn[i] - 1) * psicl(n + 2 * i + 1, n + 2 * L + 2)**j * psicl(n + 2 * i + 2, n + 2 * L + 2)**(tvn[i] - 1 - j)
                       for j in range(tvn[i])]) for i in range(L)] + [sum([(-psicl(n + 2 * L + 1, n + 2 * L + 2))**j * psicl(n + 2 * L + 2, n + 2 * L + 2)**(t1 - 1 - j) for j in range(t1)])])
            result += QQ((1, 2)) * tautclass([decstratum(graph, poly=psi)], g, n)
    return result
##
##


def generalized_hodge_chern_exp(t, a, g=None, n=None):  # The exponential e^C
    g, n = handle_g_n(g, n)
    # a.sort(key=operator.itemgetter(1,2)) ## Order a to our convention
    result = tautclass([], g, n)
    # loop on Partition of t+1 over elements of ChS, assume a is ordered
    for tvect in IntegerVectors(t + 1, len(a)):
        tvn = [tvi for tvi in tvect if tvi != 0]
        atn = [ai for ai, tvi in zip(a, tvect) if tvi != 0]
        L = len(tvn)
        Scomp = list(set(range(1, n + 1)) - set(atn[-1][1]))
        if all(atn[i][0] <= atn[i + 1][0] and atn[i][1] <= atn[i + 1][1] for i in range(L - 1)) and tvn[-1] > 1:
            vertices = [atn[0][0]] + [atn[i][0] - atn[i - 1][0] for i in range(1, L)] + [g - atn[-1][0]]
            labels = [atn[0][1] + [n + 1]] + [list(set(atn[i][1]) - set(atn[i - 1][1])) +
                                              [n + 2 * i, n + 2 * i + 1] for i in range(1, L)] + [Scomp + [n + 2 * L]]
            involution = [(n + 2 * i + 1, n + 2 * i + 2) for i in range(L)]
            graph = stgraph(vertices, labels, involution)
            psi = prod([atn[i][2]**tvn[i] * (-1)**(tvn[i] - 1) *
                        sum([binomial(tvn[i] - 1, j) * psicl(n + 2 * i + 1, n + 2 * L)**j * psicl(n + 2 * i + 2, n + 2 * L)**(tvn[i] - 1 - j) for j in range(tvn[i])]) for i in range(L - 1)] +
                       [atn[-1][2]**tvn[-1] * (-1)**(tvn[-1] - 1) *
                        sum([binomial(tvn[-1] - 1, j) * psicl(n + 2 * L - 1, n + 2 * L)**j * psicl(n + 2 * L, n + 2 * L)**(tvn[-1] - 1 - j - 1) for j in range(tvn[-1] - 1)])])
            result += multinomial(tvn) * tautclass([decstratum(graph, poly=psi)], g, n)
    return ZZ.one() / ZZ(t + 1).factorial() * result


def generalized_hodge_chern_K(t, l, a, g, n):  # The thing with K
    result = tautclass([], g, n)
    for alpha in range(t + 2):
        beta = t + 1 - alpha
        result += (divided_bernoulli_polynomial(l, beta) / ZZ(alpha).factorial()) * \
            generalized_hodge_chern_Kin(alpha, beta, a, g, n)
    return result


def generalized_hodge_chern_S(t, l, d, a, g, n):  # The thing with sigma
    result = tautclass([], g, n)
    for alpha in range(t + 1):
        beta = t + 1 - alpha
        for p in range(1, n + 1):
            result += sum(((-1)**b * d[p - 1]**(beta - b) * divided_bernoulli_polynomial(l, b) / (ZZ(alpha).factorial() * ZZ(
                beta - b).factorial())) * generalized_hodge_chern_Sin(alpha, beta, p, a, g, n) for b in range(beta + 1))
    return result


def generalized_hodge_chern_Kin(alpha, beta, a, g, n):
    if beta == 0:
        return ZZ(alpha).factorial() * generalized_hodge_chern_exp(alpha - 1, a, g, n)
    if alpha == 0:
        return kappaclass(beta - 1, g, n)
    result = tautclass([], g, n)
    # Partition of t+1 over elements of ChS, assume a is ordered
    for tvect in IntegerVectors(alpha, len(a)):
        tvn = [tvi for tvi in tvect if tvi != 0]
        atn = [ai for ai, tvi in zip(a, tvect) if tvi != 0]
        L = len(tvn)
        Scomp = list(set(range(1, n + 1)) - set(atn[-1][1]))
        if all(atn[i][0] <= atn[i + 1][0] and atn[i][1] <= atn[i + 1][1]
               for i in range(L - 1)):
            vertices = [atn[0][0]] + [atn[i][0] - atn[i - 1][0] for i in range(1, L)] + [g - atn[-1][0]]
            labels = [atn[0][1] + [n + 1]] + [list(set(atn[i][1]) - set(atn[i - 1][1])) +
                                              [n + 2 * i, n + 2 * i + 1] for i in range(1, L)] + [Scomp + [n + 2 * L]]
            involution = [(n + 2 * i + 1, n + 2 * i + 2) for i in range(L)]
            graph = stgraph(vertices, labels, involution)
            psi = prod([atn[i][2]**tvn[i] * (-1)**(tvn[i] - 1) * sum([binomial(tvn[i] - 1, j) * psicl(n + 2 * i + 1, n + 2 * L)
                       ** j * psicl(n + 2 * i + 2, n + 2 * L)**(tvn[i] - 1 - j) for j in range(tvn[i])]) for i in range(L)])
            kappa = kappacl(L, beta - 1, L + 1, g - atn[-1][0], len(Scomp) + 1)
            result += multinomial(tvn) * tautclass([decstratum(graph, poly=kappa * psi)], g, n)
    return result


def generalized_hodge_chern_Sin(alpha, beta, p, a, g, n):
    if beta == 0:
        return ZZ(alpha).factorial() * generalized_hodge_chern_exp(alpha - 1, a, g, n)
    if alpha == 0:
        return (-psiclass(p, g, n))**(beta - 1)
    result = tautclass([], g, n)
    # Partition of t+1 over elements of ChS, assume a is ordered
    for tvect in IntegerVectors(alpha, len(a)):
        tvn = [tvi for tvi in tvect if tvi != 0]
        atn = [ai for ai, tvi in zip(a, tvect) if tvi != 0]
        L = len(tvn)
        Scomp = list(set(range(1, n + 1)) - set(atn[-1][1]))
        if all(atn[i][0] <= atn[i + 1][0] and atn[i][1] <= atn[i + 1][1]
               for i in range(L - 1)) and p in Scomp:
            vertices = [atn[0][0]] + [atn[i][0] - atn[i - 1][0] for i in range(1, L)] + [g - atn[-1][0]]
            labels = [atn[0][1] + [n + 1]] + [list(set(atn[i][1]) - set(atn[i - 1][1])) +
                                              [n + 2 * i, n + 2 * i + 1] for i in range(1, L)] + [Scomp + [n + 2 * L]]
            involution = [(n + 2 * i + 1, n + 2 * i + 2) for i in range(L)]
            graph = stgraph(vertices, labels, involution)
            psi = prod([atn[i][2]**tvn[i] * (-1)**(tvn[i] - 1) * sum([binomial(tvn[i] - 1, j) * psicl(n + 2 * i + 1, n + 2 * L)**j * psicl(n + 2 *
                       i + 2, n + 2 * L)**(tvn[i] - 1 - j) for j in range(tvn[i])]) for i in range(L)] + [(-1)**(beta - 1) * psicl(p, n + 2 * L)**(beta - 1)])
            result += multinomial(tvn) * tautclass([decstratum(graph, poly=psi)], g, n)
    return result


def twist_div_to_zero(l, d, g=None, n=None):
    """
    INPUT:

    - l -- integer

    - d -- list of integers (d_1,...,d_n)

    EXAMPLES::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: twist_div_to_zero(3,[3,4],2,2)
        [[0, [1, 2], -10], [1, [1], -9], [1, [1, 2], -16]]
    """
    g, n = handle_g_n(g, n)
    seprange = ([h, [1] + s] for h in range(g + 1)
                for s in powerset(range(2, n + 1))
                if (h != 0 or s) and (h != g or len(s) <= n - 3))
    a = []
    for hS0, hS1 in seprange:
        ahS = -l * (2 * hS0 - 2 + len(hS1) + 1) - sum(d[p - 1] for p in hS1)
        if ahS != 0:
            a.append([hS0, hS1, ahS])
    return a


def generalized_lambda(deg, l, d, a, g=None, n=None):
    r"""
    Computes the Chern class c_deg of the derived pushforward of a line
    bundle \O(D) on the universal curve C_{g,n} over the space Mbar_{g,n} of stable curves.

    A divisor D on C_{g,n} up to pullbacks from Mbar_{g,n} takes the form

      D = l \tilde{K} + sum_{i=1}^n d_i \sigma_i  +  \sum_{h,S} a_{h,S} C_{h,S}

    where the numbers l, d_i and a_{h,S} are integers, \tilde{K} is the relative canonical
    class of the morphism C_{g,n} -> Mbar_{g,n}, \sigma_i is the image of the ith section
    and C_{h,S} is the boundary divisor of C_{g,n} where the moving point lies on a genus h
    component with markings given by the set S.
    For such a divisor this function computes the Chern class

      c_deg(R^\bullet \pi_* \O(D))

    of the derived pushforward R^\bullet \pi_* \O(D) using the formula given in
    [PRvZ20]_. When l=1, d=[0,...,0] and a=[] this coincides with the lambda class.

    NOTE: This formula assumes that the number n of marked points is at least 1.

    INPUT:

    deg : integer
      degree of the chern class
    l : integer
      the coefficient in front of \tilde{K}
    d : list
      vector of integers (d_1,...,d_n)
    a : list
      list of triples [h,S,a_{h,S}] where a_{h,S} is nonzero (and S is a subset of [1,...,n] )
    g : integer
      genus
    n : integer
      number of markings, for the formula to work we need n > 0


    EXAMPLES::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: l = 1; d = [0]; a = []
        sage: t = generalized_lambda(1,l,d,a,2,1);t
        Graph :      [2] [[1]] []
        Polynomial : 1/12*(kappa_1)_0 - 1/12*psi_1
        <BLANKLINE>
        Graph :      [1] [[3, 4, 1]] [(3, 4)]
        Polynomial : 1/24
        <BLANKLINE>
        Graph :      [1, 1] [[3], [1, 4]] [(3, 4)]
        Polynomial : 1/12
        sage: t.basis_vector()
        (1/2, -1/2, -1/2)

    We can verify that this computes the lambda class when l=1, d=[0,...,0] and
    a=[]. We do this for lambda_2 in \barM_3,1::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: reset_g_n(3, 1)
        doctest:...: DeprecationWarning: reset_g_n is deprecated. Please use TautologicalRing instead.
        See https://gitlab.com/modulispaces/admcycles/-/merge_requests/109 for details.
        sage: l = 1; d = [0]; a = []
        sage: s = lambdaclass(2)
        doctest:...:     DeprecationWarning: Use of lambdaclass without specifying g,n is deprecated.
        See https://gitlab.com/modulispaces/admcycles/-/merge_requests/109 for details.
        sage: t = generalized_lambda(2, l, d, a)
        sage: (s-t).simplify()
        0
    """
    g, n = handle_g_n(g, n)

    def chclass(m):
        return generalized_hodge_chern_single(m, l, d, a, g, n)

    chpoly = chern_char_to_poly(chclass, deg, g, n).degree_part(deg)
    chpoly.simplify()
    return chpoly


def DR_phi(g, d):
    r"""
    Computes the double ramification cycle DR(g,d) in the case d=[0,...,\pm
    1,...,\mp 1,...,0] as described in the paper [PRvZ20]_.

    Input:

    g : integer
      the genus
    d : vector
      the vector d

    EXAMPLES::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: DR_phi(1, [1,-1])
        Graph :      [1] [[1, 2]] []
        Polynomial : -1/12*(kappa_1)_0 + 13/12*psi_1 + 1/12*psi_2
        <BLANKLINE>
        Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
        Polynomial : -1/24
        <BLANKLINE>
        Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
        Polynomial : -1/12

    The different ways of computing the DR cycle coincide when d is of the form
    [0,...,\pm 1,...,\mp 1,...,0]::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: g = 2; d = [1,-1]
        sage: (DR_phi(g,d) - DR_cycle(g,d)).is_zero()
        True

    But not for other values of d::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: g = 1; d = [2,0]
        sage: (DR_phi(g,d)-DR_cycle(g,d)).FZsimplify()
        Graph :      [1] [[1, 2]] []
        Polynomial : -(kappa_1)_0 + 4*psi_1

    On the locus of compact type curves we have equality for all vectors d where
    the entries sum to 0.::

        sage: from admcycles import *
        sage: from admcycles.GRRcomp import *
        sage: g = 2; d = [2,-2]
        sage: (DR_cycle(g,d) - DR_phi(g,d)).basis_vector()
        (12, -4, 14, 7, -40, -10, -14, -12, 28, -4, 6, -1, 4, 0)
        sage: (DR_cycle(g,d) - DR_phi(g,d)).toTautbasis(moduli='ct')
        doctest:...: DeprecationWarning: toTautbasis is deprecated. Please use basis_vector instead.
        See https://gitlab.com/modulispaces/admcycles/-/merge_requests/109 for details.
        (0, 0, 0, 0, 0)
    """
    n = len(d)
    a = twist_div_to_zero(0, d, g, n)
    return chern_char_to_class(g, -generalized_hodge_chern(0, d, a, g, g, n))


def DR_phi1(g, d):
    r"""
    Alternative to DR_phi using chern_char_to_poly, but this takes a lot longer, kept here
    only for testing purposes.

    TESTS::

      sage: from admcycles.GRRcomp import DR_phi, DR_phi1
      sage: D1 = DR_phi(2, [1, -1])
      sage: D2 = DR_phi1(2, [1, -1])
      sage: (D1 - D2).is_empty()
      True
    """
    n = len(d)
    a = twist_div_to_zero(0, d, g, n)

    def chclass(m):
        temp = generalized_hodge_chern_single(m, 0, d, a, g, n)
        temp.simplify()
        return temp

    chpoly = chern_char_to_poly(chclass, g, g, n).degree_part(g)
    chpoly.simplify()
    return chpoly


def hodge_chern_class_matrix(char, p, q, g=None, n=None):
    g, n = handle_g_n(g, n)
    c = [chern_char_to_class(t, char, g, n) for t in range(q - p + 1, q + p)]
    return matrix([[c[i + p - 1 - j] for i in range(p)] for j in range(p)])


def generalized_brill_noether_class(r, deg, l, d, a, g=None, n=None):
    g, n = handle_g_n(g, n)
    p = r + 1
    q = g - deg + r
    char = tautclass([], g, n) + generalized_hodge_chern(l, d, a, q + p - 1, g, n)
    return hodge_chern_class_matrix(char, p, q, g, n).det()