generate.py

back to table · edit · history · where entries came from · files · download

7574 bytes, as of the version from 2026-11-02 18:57 (current). Recorded here, not run.

"""Critical L-values of higher-level newforms -- numberdb.org/T683.

Run with SageMath and the attachments roster.json and hecke-nf.json:
    sage -pip install numberdb
    sage -python generate.py                 # verify every stored entry
    sage -python generate.py --publish       # fill the existing draft with a key

The default selection is 138 embedded newforms and every critical integer,
418 values. See argument-selection.md. Exact coefficient fields, rational
isolating rectangles and coefficients through the Sturm bound identify every
embedding. PARI's list positions are never used as public identities.

Fresh PARI lfunmf/lfun computations use 160 and 220 decimal working digits;
explicit bit precision is passed to every numerical PARI call. Componentwise
agreement supports the 100 recorded digits, without a certified tail bound.
For real forms, an identically zero imaginary component is established by
exact conjugation through the Sturm bound. Forced central zeros follow from
an exact Fricke calculation. No derivative replaces a zero.
"""
import sys
import json
from pathlib import Path
from functools import lru_cache
import numberdb.sage as numberdb
from sage.libs.pari import pari
from sage.rings.rational_field import QQ
from sage.rings.integer_ring import ZZ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.qqbar import QQbar, _init_qqbar
from sage.rings.complex_mpfr import ComplexField
from sage.rings.complex_interval_field import ComplexIntervalField

# Sage's named import does not run this initializer, which exact complex
# conjugation needs. Initialize it explicitly without importing sage.all.
_init_qqbar()

BASE = Path(__file__).resolve().parent
ROSTER = json.loads((BASE / 'roster.json').read_text())
BY_FORM = {e['label']: (r, e) for r in ROSTER for e in r['embeddings']}
WORKING_DIGITS = (160, 220)
P = PolynomialRing(QQ, 'x')


def rational_decimal(text):
    from fractions import Fraction
    x = Fraction(text)
    return QQ(x.numerator) / QQ(x.denominator)


def coeff_poly(value):
    return P(str(pari.lift(value)).replace('y', 'x'))


@lru_cache(maxsize=None)
def algebra(label):
    row = next(r for r in ROSTER if r['label'] == label)
    F = P(row['field_poly'])
    cs = [P([QQ(c) for c in v]) for v in row['an_power_basis']]
    roots = [z for z, mult in F.roots(QQbar) if mult == 1]
    assert len(roots) == row['degree']
    embedded = {}
    for e in row['embeddings']:
        box = [[rational_decimal(c) for c in edge] for edge in e['root_rectangle']]
        found = [z for z in roots if box[0][0] < z.real() < box[0][1]
                 and box[1][0] < z.imag() < box[1][1]]
        assert len(found) == 1, e['label']
        z = found[0]
        is_real = all(c(z).imag() == 0 for c in cs)
        embedded[e['label']] = (z, is_real)
    zs = [z for z, real in embedded.values()]
    assert all(zs[i] != zs[j] for i in range(len(zs)) for j in range(i))
    return F, cs, embedded


def reconstruct(row):
    N, k = row['level'], row['weight']
    spec = [N, k, row['character_discriminant']]
    mf = pari.mfinit(spec, 0)
    F, cs, embedded = algebra(row['label'])
    h = P([QQ(c) for c in row['pari_generator_in_catalogue_field']])
    choices = []
    for f in pari.mfeigenbasis(mf):
        fp = pari.mfparams(f)[3]
        if str(fp) != row['pari_field_poly']:
            continue
        cf = list(pari.mfcoefs(f, row['sturm_bound']))
        if all(coeff_poly(a)(h) % F == c for a, c in zip(cf, cs)):
            choices.append(f)
    assert len(choices) == 1, ('exact reconstruction', row['label'])
    f = choices[0]
    if row['fricke_eigenvalue'] is not None:
        A = pari.mfatkininit(mf, N, precision=128)
        assert A[2] == 1
        af = pari.mfatkin(A, f)
        ac = list(pari.mfcoefs(af, row['sturm_bound']))
        fc = list(pari.mfcoefs(f, row['sturm_bound']))
        assert all(a == row['fricke_eigenvalue'] * b for a, b in zip(ac, fc))
    return mf, f, h, embedded


@lru_cache(maxsize=None)
def orbit_values(label, working):
    row = next(r for r in ROSTER if r['label'] == label)
    bits = numberdb.bits(working, losing=32)
    C = ComplexField(bits)
    pari.set_real_precision(working + 12)
    mf, f, h, embedded = reconstruct(row)
    d, k = row['degree'], row['weight']
    raw = pari.lfunmf(mf, f, precision=bits)
    Ls = [raw] if d == 1 else list(raw)
    fp = pari.mfparams(f)[3]
    images = [C(0)] if d == 1 else [C(z) for z in pari.mfembed(f, pari.Mod('y', fp), precision=bits)]
    result = {}
    used = set()
    for full, (root, is_real) in embedded.items():
        target = C(h(root))
        matches = [i for i, z in enumerate(images) if abs(z-target) < C(10).real()**(-working+20)]
        assert len(matches) == 1, (full, 'numerical root matching')
        j = matches[0]
        assert j not in used
        used.add(j)
        L = Ls[j]
        eps = C(L[5])
        if row['fricke_eigenvalue'] is not None:
            exact_eps = (-1)**(k//2) * row['fricke_eigenvalue']
            assert eps == exact_eps
        else:
            exact_eps = None
        fe_bits = int(pari.lfuncheckfeq(L, precision=bits))
        assert fe_bits < -3 * (working-20), (full, fe_bits)
        vals = {}
        for s in range(1, k):
            raw_value = C(pari.lfun(L, s, precision=bits))
            if exact_eps == -1 and 2*s == k:
                assert abs(raw_value) < C(10).real()**(-working+20)
                vals[s] = ZZ(0)
            else:
                if is_real:
                    assert abs(raw_value.imag()) < C(10).real()**(-working+20)
                    raw_value = C(raw_value.real(), 0)
                vals[s] = raw_value
        result[full] = {'values': vals, 'root_number': eps, 'is_real': is_real,
                        'functional_equation_bits': fe_bits}
    assert len(used) == d
    return result


class CriticalNewformLValues(numberdb.Generator):
    table = 'T683'
    parameters = ('form', 's')
    type = 'C'
    digits = 100
    rigour = 'heuristic (agreement-checked)'
    files = ('generate.py', 'roster.json', 'hecke-nf.json', 'checks.py',
             'check.py', 'argument-selection.md', 'numerical-evidence.json',
             'source-provenance.json', 'cc-records.json')

    def enumerate(self):
        for row in ROSTER:
            for e in row['embeddings']:
                for s in range(1, row['weight']):
                    yield {'form': e['label'], 's': s}

    def value(self, params, digits):
        assert digits <= self.digits, 'Reassess working precision before extending digits.'
        full, s = params['form'], int(params['s'])
        row, e = BY_FORM[full]
        assert 1 <= s < row['weight']
        values = {w: orbit_values(row['label'], w)[full]['values'][s] for w in WORKING_DIGITS}
        if all(v == 0 for v in values.values()):
            return ZZ(0)
        real = numberdb.agreeing(lambda w: values[w].real(), at=WORKING_DIGITS)
        if orbit_values(row['label'], WORKING_DIGITS[0])[full]['is_real']:
            return real
        imag = numberdb.agreeing(lambda w: values[w].imag(), at=WORKING_DIGITS)
        return ComplexIntervalField(numberdb.bits(max(WORKING_DIGITS)))(real, imag)


if __name__ == '__main__':
    if '--publish' in sys.argv:
        print(CriticalNewformLValues().publish(message='Compute critical L-values at two precisions with exact embedding identities'))
    else:
        report = CriticalNewformLValues().verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)