checks.py

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

5680 bytes, as of the version from 2026-11-02 17:55 (current). Recorded here, not run.

"""Independent checks for numberdb.org/T682; run: sage -python checks.py.

This is build evidence, not the campaign's separate independent critique.
The third moment uses LS section 5.1's cotangent sum, grouped by residue
class and evaluated with Hurwitz zeta, independent of the polylog root-pair
implementation in generate.py. The second moment uses exact Bernoulli root
pairs; polynomial identities use exact rational polynomials. Direct singular
quadrature uses mpmath and supplies no certification.
"""
import json
import math
import sys
from pathlib import Path

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.complex_arb import ComplexBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
import mpmath as mp
from generate import GeometricSumMahler, moment


def independent(n, k, precision=560):
    R = RealBallField(precision)
    if k == 2:
        # pi^2/2 sum_{i,j} B_2({alpha_i-alpha_j}).
        angles = [QQ(j)/n for j in range(1,n)]
        total = QQ(0)
        for a in angles:
            for b in angles:
                q = a-b
                q -= q.floor()
                total += q*q-q+QQ(1)/6
        assert total/2 == QQ(n-1)/(6*n)
        value = R.pi()**2*total/2
    else:
        C = ComplexBallField(precision)
        total = R(0)
        for r in range(1,n):
            if 2*r == n:
                continue  # cot(pi/2)=0 exactly.
            angle = R.pi()*QQ(r)/n
            hurwitz = C(2).zeta(C(QQ(r)/n))
            assert hurwitz.imag().contains_zero()
            total += angle.cos()/angle.sin()*hurwitz.real()
        value = QQ(9*n-6-3*n**3)/(2*n**2)*R(3).zeta() + 3*R.pi()*total/(2*n**2)
    assert value.is_finite() and not value.contains_zero()
    assert R(value.rad()) < abs(R(value.mid()))*R(10)**(-140)
    return value


def direct_integral(n, k, dps=45):
    # On t=(j+u)/n, |sin(n*pi*t)|=sin(pi*u), 0<u<1.
    # Evaluate each half using distance from its nearer endpoint, which avoids
    # cancellation of 1-u in tanh-sinh nodes extremely close to a root.
    with mp.workdps(dps):
        total = mp.mpf(0)
        for j in range(n):
            def f_left(v):
                if v == 0:
                    return mp.log(n)**k if j == 0 else mp.mpf(0)
                den_t = (mp.mpf(j)+v)/n
                if den_t > mp.mpf(1)/2:
                    den_t = (mp.mpf(n-j)-v)/n
                logabs = mp.log(mp.sinpi(v))-mp.log(mp.sinpi(den_t))
                return logabs**k
            def f_right(v):
                if v == 0:
                    return mp.log(n)**k if j == n-1 else mp.mpf(0)
                den_t = (mp.mpf(j+1)-v)/n
                if den_t > mp.mpf(1)/2:
                    den_t = (mp.mpf(n-j-1)+v)/n
                logabs = mp.log(mp.sinpi(v))-mp.log(mp.sinpi(den_t))
                return logabs**k
            total += mp.quad(f_left,[0,mp.mpf(1)/2])+mp.quad(f_right,[0,mp.mpf(1)/2])
        return total/n


def run():
    # Known controls before trusting either numerical method.
    with mp.workdps(45):
        for k, target in [(2,mp.pi**2/12),(3,-3*mp.zeta(3)/2)]:
            computed = direct_integral(2,k)
            assert abs(computed-target) < mp.mpf('1e-35')
            print('quadrature control n=2 k=%d passed' % k, flush=True)
    for k in (2,3):
        R=RealBallField(560)
        target=R.pi()**2/12 if k==2 else -QQ(3)/2*R(3).zeta()
        assert independent(2,k).overlaps(target)
        assert moment(2,k).overlaps(target)
    P=PolynomialRing(QQ,'x');x=P.gen()
    polynomial_checks=[]
    for n in range(2,51):
        s=sum(x**j for j in range(n))
        product=P(1)
        for divisor in range(2,n+1):
            if n % divisor==0:
                product *= P.cyclotomic_polynomial(divisor)
        assert product == s
        assert (x-1)*s == x**n-1
        if n%2:
            q=(n-1)//2
            assert s == sum(x**(q+j) for j in range(-q,q+1))
        polynomial_checks.append(n)
    results=[]
    for params in GeometricSumMahler().enumerate():
        n,k=params['n'],params['k']
        value=moment(n,k); high=moment(n,k,130); other=independent(n,k)
        assert value.overlaps(high) and value.overlaps(other)
        assert value > 0 if k==2 else value < 0
        quadrature=direct_integral(n,k)
        with mp.workdps(45):
            error=abs(quadrature-mp.mpf(str(value.mid())))
            assert error < mp.mpf('1e-34'),(n,k,str(error))
            qtext=mp.nstr(quadrature,40); etext=mp.nstr(error,5)
        results.append({'n':n,'k':k,'ball':str(value),'accuracy_bits':int(value.accuracy()),
                        'radius':str(value.rad()),'independent_ball':str(other),
                        'quadrature':qtext,'quadrature_absolute_difference':etext})
        if k==3: print('checked length',n,flush=True)
    evidence={'table':'T682','stage':'build','skill_manifest':
        'agents/runs/20261102T172655Z-build-skills/manifest.json',
        'skill_sha256':'a2b28d9953264d7b0ca62663785c2952190c5be6d6d3597fd695d9a54d606278',
        'entries':len(results),'working_bits':numberdb.bits(100,losing=128),
        'minimum_accuracy_bits':min(r['accuracy_bits'] for r in results),
        'independent_method':'Bernoulli root pairs for k=2; LS section 5.1 cotangent/Hurwitz zeta for k=3 at 560 bits',
        'quadrature_dps':45,'quadrature_tolerance':'1e-34','controls':'n=2 gamma-ratio moments, both k=2 and k=3',
        'polynomial_identities_verified_for':polynomial_checks,'results':results}
    print('EVIDENCE_JSON='+json.dumps(evidence),flush=True)


if __name__ == '__main__':
    run()