GRRcomp.py

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 [PRvZ] 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 [PRvZ]:

- ``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 [PRvZ] 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 [JPPZ]. 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 [Pagani-Ricolfi-van Zelm].

  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^1 )_0
    <BLANKLINE>
    Graph :      [1] [[1, 2]] []
    Polynomial : (-13/12)*psi_1^1
    <BLANKLINE>
    Graph :      [1] [[1, 2]] []
    Polynomial : (-1/12)*psi_2^1
    <BLANKLINE>
    Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
    Polynomial : (-11/12)*
    <BLANKLINE>
    Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
    Polynomial : 1/24*
  """
  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
  [Pagani-Ricolfi-van Zelm].

  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^1 )_0
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : (-13/12)*psi_1^1
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : (-1/12)*psi_2^1
      <BLANKLINE>
      Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
      Polynomial : (-11/12)*
      <BLANKLINE>
      Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
      Polynomial : 1/24*

  If you do not want to manually enter ``g`` and ``n`` each time these can be
  set globally::

      sage: reset_g_n(1, 2)
      sage: generalized_hodge_chern_single(1, l, d, a)
      Graph :      [1] [[1, 2]] []
      Polynomial : 1/12*(kappa_1^1 )_0
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : (-13/12)*psi_1^1
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : (-1/12)*psi_2^1
      <BLANKLINE>
      Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
      Polynomial : (-11/12)*
      <BLANKLINE>
      Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
      Polynomial : 1/24*
  """
  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)
  return result.simplify(g,n,t)




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, 1] [[1, 2, 4], [5]] [(4, 5)]
    Polynomial : 1/12*
    <BLANKLINE>
    Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
    Polynomial : 1/24*
  """
  g, n = handle_g_n(g, n)
  a.sort(key=operator.itemgetter(1, 2))  # Order a to our convention
  result = tautclass([])
  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([])

  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)]) ## 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)]) for ahS in seprange) ## The separating graphs

  result = divided_bernoulli(t+1)*(QQ((1,2)) * tautclass1 + tautclass2 )
  return   (result+0*fundclass(g,n)).simplify(g,n,t)


##
## The graphs in the sigma part
def todd_inv_exp_graphs(a,t1,tvect,g,n):
  #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 = 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)])

        # 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)])

    #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)])
  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([])
  # 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)])
  return ZZ.one()/ZZ(t+1).factorial()*result


def generalized_hodge_chern_K(t,l,a,g,n):#The thing with K
  result=tautclass([])
  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([])
  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([])
  # 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)])
  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([])
  # 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)])
  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
  [Pagani-Ricolfi-van Zelm]. 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^1 )_0
      <BLANKLINE>
      Graph :      [2] [[1]] []
      Polynomial : (-1/12)*psi_1^1
      <BLANKLINE>
      Graph :      [1, 1] [[3], [1, 4]] [(3, 4)]
      Polynomial : 1/12*
      <BLANKLINE>
      Graph :      [1] [[3, 4, 1]] [(3, 4)]
      Polynomial : 1/24*
      sage: t.toTautbasis()
      (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)
      sage: l = 1; d = [0]; a = []
      sage: s = lambdaclass(2)
      sage: t = generalized_lambda(2, l, d, a)
      sage: (s-t).simplify()
      <BLANKLINE>
  """
  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)
  chpoly.simplify(g,n,deg)
  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 [Pagani-Ricolfi-van Zelm].

  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^1 )_0
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : 13/12*psi_1^1
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : 1/12*psi_2^1
      <BLANKLINE>
      Graph :      [0, 1] [[1, 2, 4], [5]] [(4, 5)]
      Polynomial : (-1/12)*
      <BLANKLINE>
      Graph :      [0] [[4, 5, 1, 2]] [(4, 5)]
      Polynomial : (-1/24)*

  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 : (-1)*(kappa_1^1 )_0
      <BLANKLINE>
      Graph :      [1] [[1, 2]] []
      Polynomial : 4*psi_1^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)).toTautbasis()
      (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')
      (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).toTautvect().is_zero()
    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)
  chpoly.simplify(g,n,g)
  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([]) + generalized_hodge_chern(l,d,a,q+p-1,g,n)
  return hodge_chern_class_matrix(char,p,q,g,n).det()