CoCalc -- evaluation.py

Introduction lectures / admcycles / DR / evaluation.py
¹⁸¹ views
unlisted
ubuntu2004
1
r"""
2
Evaluation of integrals of cohomology classes against the fundamental class
3
"""
4

5
from sage.misc.cachefunc import cached_function
6
from sage.rings.integer_ring import ZZ
7
from sage.combinat.subset import Subsets
8
from sage.arith.misc import bernoulli, multinomial
9

10

11
from .moduli import MODULI_ST, MODULI_CT, MODULI_RT, MODULI_SM, dim_form
12
from .utils import setparts_with_auts
13
from .graph import single_stratum
14

15

16
def socle_evaluation(num, g, markings=(), moduli_type=MODULI_ST):
17
    r"""
18
    EXAMPLES::
19

20
        sage: from admcycles.DR.graph import num_strata
21
        sage: from admcycles.DR.moduli import MODULI_ST
22
        sage: from admcycles.DR.evaluation import socle_evaluation
23
        sage: g = 2
24
        sage: markings = (1,)
25
        sage: for i in range(num_strata(g, 3*g-3+len(markings), (1,))):
26
        ....:     print(i, socle_evaluation(i, g, markings, MODULI_ST))
27
        0 1/1152
28
        1 13/1920
29
        2 53/5760
30
        3 259/5760
31
        4 29/128
32
        5 1/384
33
        6 101/5760
34
        7 169/1920
35
        8 29/5760
36
        9 139/5760
37
        10 29/5760
38
        11 1/1152
39
        12 1/576
40
        ...
41
        76 2
42
        77 2
43
        78 1
44
        79 1
45
        80 1
46
        81 1
47
        82 1
48
        83 1
49
        84 1
50
        85 1
51
        86 1
52
        87 1
53
        88 1
54
        89 1
55
        90 1
56
        91 1
57
    """
58
    answer = 1
59
    G = single_stratum(num, g, dim_form(
60
        g, len(markings), moduli_type), markings, moduli_type)
61
    for i in range(1, G.M.nrows()):
62
        g0 = G.M[i, 0][0]
63
        psilist = []
64
        for j in range(1, G.M.ncols()):
65
            if G.M[i, j][0] > 0:
66
                psilist.append(G.M[i, j][1])
67
                if G.M[i, j][0] == 2:
68
                    psilist.append(G.M[i, j][2])
69
        n0 = len(psilist)
70
        dim0 = dim_form(g0, n0, moduli_type)
71
        kappalist = []
72
        for j in range(1, dim0 + 1):
73
            for k in range(G.M[i, 0][j]):
74
                kappalist.append(j)
75
        if sum(psilist) + sum(kappalist) != dim0:
76
            raise ValueError("wrong dimension")
77
        answer *= socle_formula(g0, psilist, kappalist, moduli_type)
78
    return answer
79

80

81
def socle_formula(g, psilist, kappalist, moduli_type=MODULI_ST):
82
    r"""
83
    Return the integral of a product of kappa and psi classes.
84

85
    EXAMPLES::
86

87
        sage: from admcycles.DR.evaluation import socle_formula
88
        sage: from admcycles.DR.moduli import MODULI_SM, MODULI_RT, MODULI_CT, MODULI_ST
89

90
        sage: socle_formula(4, [], [2], MODULI_SM)
91
        1/3225600
92
        sage: socle_formula(4, [], [1, 1], MODULI_SM)
93
        1/302400
94
        sage: socle_formula(3, [0], [2], MODULI_SM)
95
        1/120960
96
        sage: socle_formula(3, [0], [1, 1], MODULI_SM)
97
        1/13440
98
        sage: socle_formula(3, [1], [1], MODULI_SM)
99
        1/24192
100
        sage: socle_formula(3, [2], [], MODULI_SM)
101
        1/120960
102

103
        sage: socle_formula(3, [0], [2], MODULI_RT)
104
        1/120960
105
        sage: socle_formula(3, [0], [1, 1], MODULI_RT)
106
        1/13440
107
        sage: socle_formula(3, [1], [1], MODULI_RT)
108
        1/24192
109
        sage: socle_formula(3, [2], [], MODULI_RT)
110
        1/120960
111

112
        sage: socle_formula(4, [], [5], MODULI_CT)
113
        127/154828800
114
        sage: socle_formula(4, [], [1, 4], MODULI_CT)
115
        127/7741440
116
        sage: socle_formula(4, [], [2, 3], MODULI_CT)
117
        2159/77414400
118
        sage: socle_formula(4, [], [1, 1, 1, 1, 1], MODULI_CT)
119
        3171571/77414400
120

121
        sage: socle_formula(3, [], [6], MODULI_ST)
122
        1/82944
123
        sage: socle_formula(3, [], [1, 5], MODULI_ST)
124
        1/5760
125
        sage: socle_formula(3, [], [2, 4], MODULI_ST)
126
        971/2903040
127
        sage: socle_formula(3, [], [1, 1, 4], MODULI_ST)
128
        2173/967680
129
        sage: socle_formula(3, [], [1, 1, 1, 1, 1, 1], MODULI_ST)
130
        176557/107520
131
    """
132
    degree = sum(psilist) + sum(kappalist)
133
    n = len(psilist)
134
    if moduli_type == MODULI_SM:
135
        if degree != g - 1 + (g == 0) - (n == 0):
136
            raise ValueError("invalid degree")
137
    elif moduli_type == MODULI_RT:
138
        if degree != g - 2 + n - (g == 0):
139
            raise ValueError("invalid degree")
140
    elif moduli_type == MODULI_CT:
141
        if degree != 2 * g - 3 + n:
142
            raise ValueError("invalid degree")
143
    elif moduli_type == MODULI_ST:
144
        if degree != 3 * g - 3 + n:
145
            raise ValueError("invalid degree")
146

147
    if moduli_type == MODULI_CT or g == 0:
148
        return CTconst(g) * CTsum(psilist, kappalist)
149
    if moduli_type <= MODULI_SM or moduli_type == MODULI_RT:
150
        return RTconst(g) * RTsum(g, psilist, kappalist)
151
    if moduli_type == MODULI_ST:
152
        return STsum(psilist, kappalist)
153

154

155
def multi2(g, sigma):
156
    r"""
157
    EXAMPLES::
158

159
        sage: from admcycles.DR.evaluation import multi2
160
        sage: multi2(3, [3, 0])
161
        1
162
        sage: multi2(3, [2])
163
        1
164
        sage: multi2(3, [2, 2, 0])
165
        10
166
    """
167
    sigma.sort()
168
    if sigma[0] == 0:
169
        total = ZZ.zero()
170
        for i in range(len(sigma) - 1):
171
            sigmacopy = sigma[1:]
172
            if sigmacopy[i] > 0:
173
                sigmacopy[i] -= 1
174
                total += multi2(g, sigmacopy)
175
        return total
176
    term = ZZ(2 * g - 3 + len(sigma)).factorial()
177
    term *= ZZ(2 * g - 1).multifactorial(2)
178
    term /= ZZ(2 * g - 1).factorial()
179
    for i in sigma:
180
        term /= ZZ(2 * i - 1).multifactorial(2)
181
    return term
182

183

184
def STsum(psilist, kappalist):
185
    kappalist.sort()
186
    total = 0
187
    for i0, i1 in setparts_with_auts(kappalist):
188
        total += (-1)**(len(i0)) * i1 * \
189
            STrecur([1 + sum(j) for j in i0] + psilist)
190
    return total * (-1)**(len(kappalist))
191

192

193
def RTsum(g, psilist, kappalist):
194
    kappalist.sort()
195
    total = ZZ.zero()
196
    for i0, i1 in setparts_with_auts(kappalist):
197
        total += (-1) ** len(i0) * i1 * \
198
            multi2(g, [1 + sum(j) for j in i0] + psilist)
199
    return total * (-1)**(len(kappalist))
200

201

202
def CTsum(psilist, kappalist):
203
    kappalist.sort()
204
    total = ZZ.zero()
205
    for i0, i1 in setparts_with_auts(kappalist):
206
        total += (-1)**len(i0) * i1 * \
207
            multinomial([1 + sum(j) for j in i0] + psilist)
208
    return total * (-1)**len(kappalist)
209

210

211
def CTconst(g):
212
    """
213
    Sequence of rational numbers related to Bernoulli numbers.
214

215
    INPUT: an integer g
216

217
    OUTPUT: a rational
218

219
    EXAMPLES::
220

221
        sage: from admcycles.DR.evaluation import CTconst
222
        sage: [CTconst(g) for g in range(12)]
223
        [1,
224
         1/24,
225
         7/5760,
226
         31/967680,
227
         127/154828800,
228
         73/3503554560,
229
         1414477/2678117105664000,
230
         8191/612141052723200,
231
         16931177/49950709902213120000,
232
         5749691557/669659197233029971968000,
233
         91546277357/420928638260761696665600000,
234
         3324754717/603513268363481705349120000]
235
    """
236
    gg = 2 * ZZ(g)
237
    power = ZZ(2)**(gg - 1)
238
    return ((power - 1) * bernoulli(gg)).abs() / (power * ZZ(gg).factorial())
239

240

241
@cached_function
242
def RTconst(g):
243
    r"""
244
    Universal constant computing the intersection number
245

246
    \int_{Mbar_{g,n}} \psi_1^{g-1} \lambda_g \lambda_{g-1}.
247

248
    EXAMPLES::
249

250
        sage: from admcycles.DR.evaluation import RTconst
251
        sage: [RTconst(g) for g in range(12)]
252
        [1,
253
         1/24,
254
         1/2880,
255
         1/120960,
256
         1/3225600,
257
         1/63866880,
258
         691/697426329600,
259
         1/13284311040,
260
         3617/541999890432000,
261
         43867/64877386884710400,
262
         174611/2265559542005760000,
263
         77683/7958292791191142400]
264

265
    TESTS::
266

267
        sage: from admcycles import TautologicalRing
268
        sage: R = TautologicalRing(1,1)
269
        sage: (R.psi(1)^0*R.lambdaclass(1)*R.lambdaclass(0)).evaluate()
270
        1/24
271
        sage: R = TautologicalRing(2,1)
272
        sage: (R.psi(1)^1*R.lambdaclass(2)*R.lambdaclass(1)).evaluate()
273
        1/2880
274
        sage: R = TautologicalRing(3,1)
275
        sage: (R.psi(1)^2*R.lambdaclass(3)*R.lambdaclass(2)).evaluate()
276
        1/120960
277
    """
278
    if g == 0:
279
        return ZZ(1)
280
    return (bernoulli(2 * ZZ(g))).abs() / (2**(2 * g - 1) * ZZ(2 * g - 1).multifactorial(2) * 2 * g)
281

282

283
def STrecur(psi):
284
    """
285
    Integral of psi classes.
286

287
    INPUT: a sorted tuple of nonegative integers
288

289
    OUTPUT: a rational
290

291
    EXAMPLES::
292

293
        sage: from admcycles.DR.evaluation import STrecur
294
        sage: STrecur((0, 0, 0))
295
        1
296
        sage: STrecur((1,))
297
        1/24
298
        sage: STrecur((4,))
299
        1/1152
300
        sage: STrecur((7,))
301
        1/82944
302
        sage: STrecur((2,) * 3)
303
        7/240
304
        sage: STrecur((2,) * 6)
305
        1225/144
306
        sage: STrecur((1,) * 3)
307
        1/12
308
    """
309
    return STrecur_calc(tuple(sorted(psi)))
310

311

312
@cached_function
313
def STrecur_calc(psi):
314
    if not psi:
315
        return ZZ.one()
316
    s = sum(psi)
317
    n = len(psi)
318
    if (s - n) % 3:
319
        return ZZ.zero()
320
    if psi[0] == 0:
321
        if s == 0 and n == 3:
322
            return ZZ.one()
323
        total = ZZ.zero()
324
        psi_end = list(psi[1:])
325
        for i in range(n - 1):
326
            psicopy = list(psi_end)
327
            if psicopy[i] > 0:
328
                psicopy[i] -= 1
329
                total += STrecur(psicopy)
330
        return total
331

332
    g = ZZ(s - n) // 3 + 1  # in ZZ
333
    d = ZZ(psi[-1])
334
    total = ZZ.zero()
335

336
    psicopy = [0, 0, 0, 0] + list(psi)
337
    psicopy[-1] += 1
338
    total += (2 * d + 3) / 12 * STrecur(psicopy)
339

340
    psicopy = [0, 0, 0] + list(psi)
341
    total -= (2 * g + n - 1) / 6 * STrecur(psicopy)
342

343
    for I in Subsets(list(range(n - 1))):
344
        # 3g - 3 + n
345
        nI = len(I) + 2
346
        degI = sum(psi[i] for i in I)
347
        if (degI - nI) % 3 != 0:
348
            continue
349
        gI = (degI - nI) // 3 + 1
350
        if 2 * gI - 2 + nI <= 0:
351
            continue
352
        psi3 = [0, 0] + [psi[i] for i in I]
353
        x = STrecur(psi3)
354
        if not x:
355
            raise ValueError
356
        psi1 = [0, 0] + [psi[i] for i in range(n - 1) if i not in I] + [d + 1]
357
        psi2 = list(psi1[1:])
358
        psi2[-1] = d
359
        total += ((2 * d + 3) * STrecur(psi1) -
360
                  (2 * g + n - 1) * STrecur(psi2)) * x
361
    total /= (2 * g + n - 1) * (2 * g + n - 2)
362
    return total
363

364
Product

Resources

Company
