back to table · edit · history · where entries came from · files · download
6324 bytes, as of the version from 2026-11-02 04:21 (current). Recorded here, not run.
"""Mathieu and umbral moonshine coefficients -- numberdb.org/T668.
Run with SageMath (or passagemath) and the numberdb package:
sage -pip install numberdb
sage -python generate.py # verify every stored entry
sage -python checks.py # independent published controls
The generator's publish() method fills the existing draft only when deliberately
called by its owner; this script never offers or makes a draft public.
Exact formal series in QQ[y,y^-1] implement Cheng–Duncan–Harvey,
arXiv:1204.2779v3, equations (2.18)–(2.25), (2.49)–(2.51), (2.62)–(2.64).
The antisymmetric theta basis is theta_-r - theta_r, so extracting y^r
requires a minus sign. The later module convention is DGO 1503.01472 (1.4).
Read argument-selection.md for the six cases and the noncanonical grade bound.
"""
import functools
import sys
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.laurent_polynomial_ring import LaurentPolynomialRing
L = LaurentPolynomialRing(QQ, 'y')
y = L.gen()
SYSTEMS = (('A1^24', 2), ('A2^12', 3), ('A3^8', 4),
('A4^6', 5), ('A6^4', 7), ('A12^2', 13))
MAX_EXPONENT = 30
def mul(a, b, degree):
"""Truncated series multiplication, with no lazy Sage series methods."""
out = [L.zero() for _ in range(degree + 1)]
for i, ai in enumerate(a[:degree + 1]):
if ai:
for j, bj in enumerate(b[:degree + 1 - i]):
if bj:
out[i + j] += ai * bj
return out
def linear(*terms):
return [sum((scale * a[n] for scale, a in terms), L.zero())
for n in range(len(terms[0][1]))]
def scalar_divide(a, denominator):
"""Formal division by a scalar-coefficient series with nonzero constant."""
out = []
for n in range(len(a)):
out.append((a[n] - sum((denominator[k] * out[n-k]
for k in range(1, n+1)), L.zero()))
/ QQ(denominator[0]))
return out
def jacobi_forms(degree):
"""Z^(m)=2 phi_1^(m), exact q-Laurent expansions through q^degree."""
# theta2 / (q^(1/8)*y^(1/2)); squaring restores a factor y.
t2 = [L.zero() for _ in range(degree + 1)]
for j in range(-degree - 2, degree + 3):
n = j*(j+1)//2
if n <= degree:
t2[n] += y**j
d2 = [f(QQ(1)) for f in mul(t2, t2, degree)]
a = scalar_divide([y*f for f in mul(t2, t2, degree)], d2)
# theta3 and theta4 use t=q^(1/2); theta4(t)=theta3(-t).
dt = 2*degree
t3 = [L.zero() for _ in range(dt + 1)]
for j in range(-dt - 1, dt + 2):
if j*j <= dt:
t3[j*j] += y**j
sq = mul(t3, t3, dt)
b = scalar_divide(sq, [f(QQ(1)) for f in sq])
c = [(-1)**n * f for n, f in enumerate(b)]
bc = mul(b, c, dt)
assert all(not bc[n] for n in range(1, dt+1, 2))
d = [b[2*n] + c[2*n] for n in range(degree+1)]
e = [bc[2*n] for n in range(degree+1)]
phi = {}
phi[2] = linear((4, a), (4, d))
phi[3] = linear((2, mul(a, d, degree)), (2, e))
phi[4] = linear((4, mul(a, e, degree)))
phi[5] = linear((QQ(1)/4, mul(phi[4], phi[2], degree)),
(-QQ(1)/4, mul(phi[3], phi[3], degree)))
phi[7] = linear((1, mul(phi[3], phi[5], degree)),
(-1, mul(phi[4], phi[4], degree)))
phi9 = linear((1, mul(phi[3], phi[7], degree)),
(-1, mul(phi[5], phi[5], degree)))
phi[13] = linear((1, mul(phi[5], phi9, degree)),
(-2, mul(phi[7], phi[7], degree)))
return {m: [2*f for f in series] for m, series in phi.items()}
def mu_polynomial_part(m, degree):
"""mu_m - (y+1)/(y-1), expanding in |q|<|y|<1."""
out = [L.zero() for _ in range(degree + 1)]
for k in range(1, degree + 1):
base = m*k*k
if base > degree:
break
out[base] += y**(-2*m*k) - y**(2*m*k)
for j in range(1, (degree-base)//k + 1):
out[base + k*j] += 2*(y**(-2*m*k-j)-y**(2*m*k+j))
return out
def divide_y_minus_one(f):
"""Exact Laurent polynomial division; evaluate the remainder explicitly."""
if not f:
return L.zero()
assert f(QQ(1)) == 0
lo = int(f.valuation())
hi = int(f.degree())
result = L.zero()
carry = QQ(0)
for j in range(lo, hi):
carry -= f[j]
result += carry*y**j
assert result*(y-1) == f
return result
@functools.lru_cache(maxsize=None)
def expansions(degree):
z = jacobi_forms(degree)
mu1 = mu_polynomial_part(1, degree)
finite = {}
for _, m in SYSTEMS:
chi = QQ(24)/QQ(m-1)
assert [f(QQ(1)) for f in z[m]] == [chi]+[0]*degree
polynomial = mul(mu1, z[m], degree)
mum = mu_polynomial_part(m, degree)
finite[m] = []
for n in range(degree+1):
f = z[m][n] - (chi if n == 0 else 0)
finite[m].append((y+1)*divide_y_minus_one(f)
+ polynomial[n] - chi*mum[n])
return z, finite
def coefficient(m, r, s, degree=None):
if s <= 0 or (s+r*r) % (4*m) or not 1 <= r < m:
raise ValueError('Invalid positive umbral grade')
n = (s+r*r)//(4*m)
_, finite = expansions(max(n, degree or 0))
answer = -finite[m][n][r]
assert answer.denominator() == 1
return ZZ(answer)
class UmbralCoefficients(numberdb.Generator):
table = 'T668'
files = ('generate.py', 'checks.py', 'check.py', 'published-controls.json',
'argument-selection.md', 'extract_controls.py', 'sources.json',
'skill-manifest.json', 'numerical-evidence.json')
parameters = ('X', 'r', 's')
type = 'Z'
rigour = 'exact'
digits = 100
def enumerate(self):
for name, m in SYSTEMS:
for r in range(1, m):
for s in range((-r*r-1) % (4*m)+1, 4*m*MAX_EXPONENT+1, 4*m):
yield {'X': name, 'r': r, 's': s}
def value(self, params, digits):
m = dict(SYSTEMS)[params['X']]
r, s = int(params['r']), int(params['s'])
degree = MAX_EXPONENT + max((mm-1)**2//(4*mm) for _, mm in SYSTEMS)
return coefficient(m, r, s, degree)
if __name__ == '__main__':
report = UmbralCoefficients().verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)