back to table · edit · history · where entries came from · files · download
9491 bytes, as of the version from 2026-10-24 20:04 (current). Recorded here, not run.
"""Feynman periods of primitive log-divergent phi4 graphs -- numberdb.org/T464.
Run with SageMath (PARI >= 2.17) and numberdb installed:
sage -pip install numberdb
sage -python generate.py # verify all stored entries
sage -python checks.py # compute and check without writing
The attached generator recomputes every non-FAIL period in the version-2
Panzer--Schnetz ancillary file, whose SHA-256 is pinned below. Download
Periods beside this file, or let the generator fetch the pinned public URL.
The Maple data is parsed with a restricted expression evaluator, never exec.
Rationals stay exact. PARI zetamult uses descending indices (unlike Sage's
Multizeta wrapper); polylogmult([s1,...],[xi,1,...]) handles the colored sums.
Every value is computed at 150 and 200 decimal digits, compared with the
source's 100 significant digits, then written to 100 significant digits.
Agreement is heuristic: no interval proof of these evaluations is asserted.
The nonreal root-of-unity arguments are displaced inward by 16 binary units
at each working precision to avoid a spurious domain refusal from rounding.
This is included in the two-precision agreement check; it is not an error proof.
The one-loop integral is exactly 1 and is returned as a rational.
Duplicate exact source expressions retain all census labels and carry equals.
"""
import ast
from functools import lru_cache
import hashlib
from pathlib import Path
import re
import sys
import urllib.request
import numberdb.sage as numberdb
from sage.libs.pari import pari
from sage.rings.complex_mpfr import ComplexField
from sage.rings.rational_field import QQ
SOURCE_URL = 'https://arxiv.org/src/1603.04289v2/anc/Periods'
SOURCE_SHA256 = '093d458994903b3b397f3ce3cd5cc9578221c065e5a6b3fb0947abf1b2ddce19'
WORKING_DIGITS = (150, 200)
# Fifty guard digits at the lower precision are measured in checks.py.
BASE = Path(__file__).resolve().parent
sys.setrecursionlimit(max(sys.getrecursionlimit(), 20000))
def split_fields(text):
depth, start, out = 0, 0, []
for i, char in enumerate(text):
if char in '[{(':
depth += 1
elif char in ']})':
depth -= 1
elif char == ',' and depth == 0:
out.append(text[start:i].strip())
start = i + 1
if depth < 0:
raise ValueError('Unbalanced source record')
if depth:
raise ValueError('Unbalanced source record')
out.append(text[start:].strip())
return out
@lru_cache(None)
def records():
path = BASE / 'Periods'
raw = path.read_bytes() if path.exists() else urllib.request.urlopen(SOURCE_URL, timeout=120).read()
if hashlib.sha256(raw).hexdigest() != SOURCE_SHA256:
raise ValueError('Ancillary source has changed; inspect it before proceeding')
out = []
pattern = r'Period\[(\d+),(\d+)\]:=\[(.*?)\]:'
for match in re.finditer(pattern, raw.decode(), re.S):
fields = split_fields(match[3])
if len(fields) != 7:
raise ValueError('Expected seven fields')
loops, index = int(match[1]), int(match[2])
out.append({'loops': loops, 'index': index, 'edges': fields[0],
'expression': fields[2], 'source_number': fields[3],
'ancestor': fields[5]})
if len(out) != 10119:
raise ValueError('Unexpected census size')
return out
@lru_cache(None)
def known_records():
return tuple(r for r in records() if r['source_number'] != 'FAIL')
@lru_cache(None)
def record_map():
return {(r['loops'], r['index']): r for r in known_records()}
@lru_cache(None)
def expression_tree(text):
return ast.parse(text.replace('^', '**'), mode='eval').body
class Evaluator:
"""Only integers, rational arithmetic, Pi, I, and zeta subscripts."""
def __init__(self, digits):
self.bits = numberdb.bits(digits)
self.C = ComplexField(self.bits)
self.atoms = {}
def evaluate(self, node):
if isinstance(node, ast.Constant) and type(node.value) is int:
return QQ(node.value)
if isinstance(node, ast.Name) and node.id in ('Pi', 'I'):
return self.C.pi() if node.id == 'Pi' else self.C.gen()
if isinstance(node, ast.UnaryOp) and isinstance(node.op, (ast.UAdd, ast.USub)):
v = self.evaluate(node.operand)
return v if isinstance(node.op, ast.UAdd) else -v
if isinstance(node, ast.BinOp):
a, b = self.evaluate(node.left), self.evaluate(node.right)
if isinstance(node.op, ast.Add): return a + b
if isinstance(node.op, ast.Sub): return a - b
if isinstance(node.op, ast.Mult): return a * b
if isinstance(node.op, ast.Div): return a / b
if isinstance(node.op, ast.Pow):
if b in QQ and QQ(b).denominator() == 1:
return a ** int(b)
return self.C(a) ** self.C(b)
if isinstance(node, ast.Subscript) and isinstance(node.value, ast.Name) and node.value.id == 'zeta':
key = ast.dump(node)
if key in self.atoms:
return self.atoms[key]
args = node.slice.elts if isinstance(node.slice, ast.Tuple) else [node.slice]
if isinstance(args[0], ast.List):
pair = args[0].elts
if len(pair) != 2: raise ValueError('Unknown colored zeta')
xi = self.evaluate(pair[0])
if self.C(xi).imag() != 0:
# PARI rejects a unit root rounded just outside its domain.
# Approach its continuous boundary value from inside; the
# 16-ulp displacement vanishes as working precision grows.
xi = self.C(xi) * (1 - QQ(2)**(-self.bits + 4))
if QQ(xi.real())**2 + QQ(xi.imag())**2 >= 1:
raise ValueError('Colored argument was not rounded inward')
weights = [int(self.evaluate(pair[1]))] + [int(self.evaluate(n)) for n in args[1:]]
value = self.C(pari.polylogmult(weights, [xi] + [1] * (len(weights)-1), precision=self.bits))
else:
weights = [int(self.evaluate(n)) for n in args]
if weights[0] < 2 or min(weights) < 1:
raise ValueError('Divergent MZV')
value = self.C(pari.zetamult(weights, precision=self.bits))
self.atoms[key] = value
return value
raise ValueError('Unsupported expression node: ' + ast.dump(node)[:120])
@lru_cache(None)
def evaluator(digits):
return Evaluator(digits)
@lru_cache(None)
def evaluate_expression(text, digits):
ev = evaluator(digits)
result = ev.C(ev.evaluate(expression_tree(text)))
if any(x.is_NaN() or x.is_infinity() for x in (result.real(), result.imag())):
raise ValueError('Nonfinite period')
# Imaginary cancellation is expected for conjugate root-of-unity sums.
if abs(result.imag()) > max(ev.C(1).real(), abs(result.real())) * QQ(10)**(-digits+25):
raise ValueError('Imaginary remainder exceeds the guard')
return result.real()
def decimal_rational(text):
"""The exact rational at the source's written decimal midpoint."""
if '.' not in text:
return QQ(text), QQ(0)
whole, fraction = text.split('.')
return QQ(int(whole + fraction)) / QQ(10)**len(fraction), QQ(10)**(-len(fraction))
@lru_cache(None)
def first_labels():
first = {}
for row in known_records():
first.setdefault(row['expression'], (row['loops'], row['index']))
return first
class FeynmanPhi4Periods(numberdb.Generator):
table = 'T464'
parameters = ('loops', 'index')
type = 'R'
digits = 100
rigour = 'heuristic (agreement-checked)'
def enumerate(self):
for row in known_records():
yield {'loops': row['loops'], 'index': row['index']}
def value(self, params, digits):
if digits != 100:
raise ValueError('This source is checked at 100 significant digits')
key = (int(params['loops']), int(params['index']))
row = record_map()[key]
if key == (1, 1):
if row['expression'] != '1': raise ValueError('One-loop control changed')
return {'number': QQ(1), 'comment': r'The one-loop period is exactly $1$, since $\int_0^\infty (1+x)^{-2}\,dx=1$.'}
low, high = [evaluate_expression(row['expression'], d) for d in WORKING_DIGITS]
source, unit = decimal_rational(row['source_number'])
if abs(high - source) > unit:
raise ValueError('Source disagrees at %r' % (key,))
if abs(high - low) > abs(high) * QQ(10)**(-120):
raise ValueError('Fewer than 120 stable digits at %r' % (key,))
value = numberdb.agreeing(lambda d: evaluate_expression(row['expression'], d), at=WORKING_DIGITS)
entry = {'number': value}
first = first_labels()[row['expression']]
if first != key:
entry['equals'] = 'HREF{#%d,%d}' % first
if key == (3, 1):
entry['comment'] = r'The wheel with three spokes has period $6\zeta(3)$; its completion is $K_5$.'
return entry
if __name__ == '__main__':
import os
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
report = FeynmanPhi4Periods().verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)