checks.py

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

7730 bytes, as of the version from 2026-10-24 20:04 (current). Recorded here, not run.

"""Independent checks for Feynman phi4 periods -- numberdb.org/T464.

Run: sage -python checks.py
Requires generate.py and check.py beside this file. Source retrieval is handled
by generate.py. To check stored entries, set NUMBERDB_CHECK_STORED=1 and give
NUMBERDB_API_KEY (or NUMBERDB_KEY_FROM_STDIN=1 and the key on stdin).
Output includes all measurements and a final JSON summary. Nothing is sent.
"""
import ast
from collections import Counter
from itertools import combinations
import json
from math import factorial
import os
import sys
import time
import generate as g
from check import exactness, measure, names_its_rings
from numberdb._write import to_text
from sage.rings.real_arb import RealBallField
from sage.rings.real_mpfr import RealField
from sage.rings.rational_field import QQ
from sage.graphs.graph import Graph


def controls():
    for digits in g.WORKING_DIGITS:
        ev = g.evaluator(digits)
        R, B = RealField(ev.bits), RealBallField(ev.bits)
        tol = R(10)**(-digits+3)
        # Known nonzero answers; these are not comparisons of two zeros.
        assert B(3).zeta() > 1
        assert abs(ev.C(g.pari.zetamult([2, 1], precision=ev.bits)).real() - R(B(3).zeta())) < tol
        assert abs(g.evaluate_expression('zeta[2]', digits)-R(B.pi()**2/6)) < tol
        assert abs(g.evaluate_expression('zeta[[-1,2]]', digits)+R(B.pi()**2/12)) < tol
        assert abs(8*g.evaluate_expression('zeta[[-1,2],1]', digits)-R(B(3).zeta())) < tol
        xi = (ev.C(1)+ev.C.gen()*ev.C(3).sqrt())/2
        xin = xi*(1-QQ(2)**(-ev.bits+4))
        li1 = ev.C(g.pari.polylogmult([1], [xin], precision=ev.bits))
        assert abs(li1+(1-xi).log()) < tol
        li2 = ev.C(g.pari.polylogmult([2], [xin], precision=ev.bits))
        assert abs(li2.real()-R(B.pi()**2/36)) < tol
    print('PASS: MZV convention, Euler relation, alternating sums, sixth-root controls', flush=True)


def graph_checks(rows):
    graphs = {}
    min_cut = 100
    for r in rows:
        edges = [tuple(sorted(x)) for x in ast.literal_eval(r['edges'])]
        vertices = sorted({v for edge in edges for v in edge})
        n = len(vertices)
        assert n == r['loops']+2
        assert len(edges)-n+1 == r['loops']+3
        assert set(Counter(v for e in edges for v in e).values()) == {4}
        assert all(a != b for a,b in edges)
        # In a 4-regular completion, no edge cut of size <=4 may separate
        # two sets both having at least two vertices. Enumerate each cut once.
        bits = {v:1<<i for i,v in enumerate(vertices)}
        bit_edges = [(bits[a],bits[b]) for a,b in edges]
        for mask in range(1,1<<(n-1)):
            pop = mask.bit_count()
            if pop < 2 or pop > n-2: continue
            cut = sum(bool(mask&a) != bool(mask&b) for a,b in bit_edges)
            assert cut >= 6, (r['loops'],r['index'],mask,cut)
            min_cut = min(min_cut,cut)
        graphs[(r['loops'],r['index'])] = edges
    # Identify the zig-zag graphs by their completed circulants, independently
    # of the numeric formula or the assumption that census index 1 is special.
    for ell in range(3,12):
        n=ell+2
        circulant=Graph([list(range(n)),[(v,(v+d)%n) for v in range(n) for d in (1,2)]],format='vertices_and_edges')
        actual=Graph(graphs[(ell,1)],format='list_of_edges')
        assert actual.is_isomorphic(circulant)
    assert len(set(graphs[(3,1)])) == 10
    k44=Graph([(a,b) for a in range(4) for b in range(4,8)],format='list_of_edges')
    assert Graph(graphs[(6,4)],format='list_of_edges').is_isomorphic(k44)
    print('PASS: all completed graphs are 4-regular and satisfy the nontrivial-cut condition; nine zig-zags identified',flush=True)
    return {'graph_count':len(graphs),'minimum_nontrivial_cut':min_cut,'zigzag_graphs':9}


def number_checks(entries):
    # Written decimals mean midpoint +/- one unit of the last decimal place.
    B = RealBallField(768)
    values={}
    for key,entry in entries.items():
        text=entry if isinstance(entry,str) else to_text(entry['number'],100)
        mid,unit=g.decimal_rational(text)
        assert mid>0
        values[key]=(mid,unit)
    for ell in range(3,12):
        coeff = QQ(4*factorial(2*ell-2))/QQ(factorial(ell)*factorial(ell-1))
        coeff *= 1-QQ(1-(-1)**ell)/QQ(2)**(2*ell-3)
        want=B(coeff)*B(2*ell-3).zeta()
        mid,unit=values[(ell,1)]
        assert want.is_finite() and B(want.rad())<B(QQ(10)**(-180))
        assert want>B(mid-unit) and want<B(mid+unit)
    # A different relation of periods, PS equation (5.3), tested with
    # propagated uncertainty of the written values, not a near-zero midpoint.
    def ball(key):
        m,u=values[key]
        return B(m).add_error(B(u))
    residual=ball((6,4))-12*ball((4,1))*ball((3,1))+16*ball((6,3))
    assert residual.is_finite() and B(residual.rad())<B(QQ(10)**(-90))
    assert residual.contains_zero()
    assert values[(1,1)] == (QQ(1),QQ(0))
    print('PASS: exact one-loop value, all nine zig-zag formulas using Arb, and K_4,4 period relation',flush=True)
    return {'zigzag_ball_checks':9,'one_loop_exact':True,'K44_relation':True}


def main():
    start=time.monotonic()
    controls()
    generator=g.FeynmanPhi4Periods()
    assert not names_its_rings(str(g.BASE/'generate.py'))
    rows=g.known_records()
    graph_report=graph_checks(rows)
    stored=os.environ.get('NUMBERDB_CHECK_STORED')=='1'
    if stored:
        if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
            g.numberdb.configure(api_key=sys.stdin.read().strip())
        doc=g.numberdb.table('T464')
        raw=doc['Numbers']
        if isinstance(raw,list):
            entries={(int(r['params']['loops']),int(r['params']['index'])):r['number'] for r in raw}
        else:
            entries={(int(ell),int(k)):(v['number'] if isinstance(v,dict) else v) for ell,rs in raw.items() for k,v in rs.items()}
        assert set(entries)==set(g.record_map())
    else:
        entries={(p['loops'],p['index']):generator.value(p,100) for p in generator.enumerate()}
    typed={k:(QQ(v) if isinstance(v,str) and '.' not in v and 'e' not in v.lower() else v) for k,v in entries.items()}
    assert not exactness(typed)
    written={k: (v if isinstance(v,str) else dict(v,number=to_text(v['number'],100))) for k,v in entries.items()}
    measured=measure(written)
    numeric_report=number_checks(entries)
    min_stable=10000; max_ulps=0; equals_count=0
    canonical={}
    for r in rows:
        key=r['loops'],r['index']
        if key==(1,1):continue
        low,high=[g.evaluate_expression(r['expression'],d) for d in g.WORKING_DIGITS]
        delta=abs(low-high)/abs(high)
        stable=float(-delta.log10()) if delta else 150.0
        min_stable=min(min_stable,stable)
        source,unit=g.decimal_rational(r['source_number'])
        ulps=float(abs(high-source)/unit); assert ulps<=1
        max_ulps=max(max_ulps,ulps)
        entry=entries[key]
        text=entry if isinstance(entry,str) else to_text(entry['number'],100)
        mid,out_unit=g.decimal_rational(text)
        assert abs(high-mid)<=out_unit
        if r['expression'] in canonical:
            assert text==canonical[r['expression']]
            equals_count+=1
        else:
            canonical[r['expression']]=text
    result={'stored':stored,'entries':len(entries),'measurement':measured,
            'graph_checks':graph_report,'independent_numbers':numeric_report,
            'minimum_stable_significant_digits':min_stable,'maximum_source_ulps':max_ulps,
            'duplicate_expression_links':equals_count,'seconds':round(time.monotonic()-start,3)}
    print('REPORT '+json.dumps(result),flush=True)
    return result


if __name__=='__main__':
    main()