generate.py

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

4638 bytes, as of the version from 2026-10-29 11:13 (current). Recorded here, not run.

"""Witten zeta values of simple Lie algebras -- numberdb.org/T607.

Run with SageMath or passagemath and the numberdb client:
    sage -pip install numberdb
    sage -python generate.py                 # verify all stored entries
    sage -python generate.py --values        # compute without website writes

This generator never creates, offers or publishes a table. The run's fill
script sends its entries to the existing private draft. For draft verification
set NUMBERDB_API_KEY in the environment.

Exact rational coefficients of pi**(2*m*N) are computed by partial fractions
of the Weyl-chamber lattice sum. See derivation.md for the proof and checks.py
for independent published formulas and direct positive-sum checks. No source
coefficient table is used to generate the values. The 96-bit precision guard
is measured in numerical-evidence.json, not inferred from decimal agreement.
"""
import sys
from functools import lru_cache
from math import comb, factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

TYPES = ('A1', 'A2', 'B2', 'G2')
MAX_M = 4
WORKING_GUARD = 96
# N positive coroots, M product at rho, |W|, poles of the a/b rational function,
# and the product of its leading coefficients. The b factor is omitted.
ROOT_DATA = {
    'A2': (3, 2, 6, ('0', '-1'), 1),
    'B2': (4, 6, 8, ('0', '-1', '-2'), 1),
    'G2': (6, 120, 12, ('0', '-1', '-2', '-3', '-3/2'), 2),
}

@lru_cache(None)
def bernoulli(n):
    """Exact Bernoulli numbers, with B_1 = -1/2, from their defining recurrence."""
    if n == 0:
        return QQ(1)
    return -sum(QQ(comb(n + 1, j)) * bernoulli(j) for j in range(n)) / QQ(n + 1)

@lru_cache(None)
def zeta_coefficient(k):
    """zeta(k)/pi**k for even k >= 0 (including zeta(0)=-1/2)."""
    assert k >= 0 and k % 2 == 0
    if k == 0:
        return -QQ(1)/2
    return (-1)**(k//2+1) * bernoulli(k) * QQ(2)**(k-1) / QQ(factorial(k))

@lru_cache(None)
def partial_fractions(algebra, s):
    """Return c_(r,k) with f(x) = sum c_(r,k)/(x-r)**k, exactly."""
    _, _, _, strings, leading = ROOT_DATA[algebra]
    poles = tuple(QQ(x) for x in strings)
    result = []
    for r in poles:
        series = [QQ(leading)**(-s)] + [QQ(0)]*(s-1)
        for t in poles:
            if r == t:
                continue
            factor = [QQ((-1)**j * comb(s+j-1, j)) * (r-t)**(-s-j)
                      for j in range(s)]
            series = [sum(series[i]*factor[j-i] for i in range(j+1))
                      for j in range(s)]
        result.extend((r,k,series[s-k]) for k in range(1,s+1))
    return tuple(result)

@lru_cache(None)
def coefficient(algebra, m):
    """Return exact (c, exponent) for Z_g(2m) = c*pi**exponent."""
    s = 2*int(m)
    assert s >= 2
    if algebra == 'A1':
        return zeta_coefficient(s), s
    n, normalization, weyl_order, strings, _ = ROOT_DATA[algebra]
    w = n*s
    poles = tuple(QQ(x) for x in strings)
    total = QQ(0)
    for r,k,c in partial_fractions(algebra,s):
        if k % 2 == 0:
            parity = QQ(1)
            if r.denominator() == 2:
                even_b = QQ(2)**(-(w-k))
                parity = even_b + (QQ(2)**k-1)*(1-even_b)
            else:
                assert r.denominator() == 1
            total += 4*c*zeta_coefficient(k)*zeta_coefficient(w-k)*parity
        # Remove all other forbidden integer a=q*b. q*b is an integer
        # exactly when the reduced denominator of q divides b.
        for q in poles:
            if q != r:
                total -= 2*c*(q-r)**(-k)*QQ(q.denominator())**(-w)*zeta_coefficient(w)
    return total * QQ(normalization)**s / QQ(weyl_order), w

class WittenZeta(numberdb.Generator):
    table = 'T607'
    parameters = ('algebra', 'm')
    type = 'R'
    digits = 100
    rigour = 'proven'

    def enumerate(self):
        for algebra in TYPES:
            for m in range(1, MAX_M+1):
                yield {'algebra':algebra,'m':str(m)}

    def value(self, params, digits):
        c,w = coefficient(params['algebra'],int(params['m']))
        field = RealBallField(numberdb.bits(digits,losing=WORKING_GUARD))
        value = field(c)*field.pi()**w
        assert value.is_finite() and value > 1
        assert value.accuracy() >= numberdb.bits(digits)
        return value

if __name__ == '__main__':
    g = WittenZeta()
    if '--values' in sys.argv:
        for params in g.enumerate():
            print(params, coefficient(params['algebra'],int(params['m'])),g.value(params,100))
    else:
        report = g.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)