generate.py

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

4616 bytes, as of the version from 2026-10-31 21:13 (current). Recorded here, not run.

"""Squared Nyman–Beurling approximation distances -- numberdb.org/T641.

Run with SageMath (or a compatible passagemath installation):
    sage -pip install numberdb
    sage -python generate.py                  # verify the stored draft/table
    sage -python generate.py --fill           # write entries, never publish it
Authentication for a private draft uses NUMBERDB_API_KEY. The campaign wrapper
may supply the key on stdin with NUMBERDB_KEY_FROM_STDIN=1.

The sole parameter is the positive integer truncation index. By default all
N=1,...,30 are computed in increasing order. The unanchored basis is
rho_n(x)={1/(n*x)} in L^2((0,infinity),dx). The target is 1_(0,1).
All arithmetic is Arb real ball arithmetic, including the LDL^T elimination.
Every pivot must be strictly positive. Therefore elimination encloses the
exact minimization problem and never divides by an interval containing zero.
The 256 extra working bits are measured in numerical-evidence.json; they are
not a replacement for the propagated enclosure. No stored rounded companion
values are used. See derivation.md and argument-selection.md.
"""
import os
import sys
from functools import lru_cache
from math import gcd

import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField

WORKING_GUARD = 256
MAX_N = 30


def gram_system(limit, precision):
    """Fresh finite cotangent sums, Gram entries, and target inner products."""
    R = RealBallField(precision)
    C = (2*R.pi()).log() - R.euler_constant()

    @lru_cache(None)
    def vasyunin(p, q):
        return sum((R(QQ((m*p) % q)/QQ(q)) *
                    (R.pi()*R(QQ(m)/QQ(q))).cot()
                    for m in range(1, q)), R(0))

    @lru_cache(None)
    def autocorrelation(p, q):
        x = R(QQ(p)/QQ(q))
        return (((1-x)*x.log() + (1+x)*C)/2
                - R.pi()*(vasyunin(p,q)+vasyunin(q,p))/(2*q))

    G = [[R(0) for k in range(limit)] for j in range(limit)]
    for j in range(1, limit+1):
        for k in range(1, j+1):
            g = gcd(j,k)
            G[j-1][k-1] = G[k-1][j-1] = autocorrelation(j//g,k//g)/j
    b = [(R(n).log()+1-R.euler_constant())/n for n in range(1,limit+1)]
    return R, G, b


def ldl_distances(G, b):
    """All nested minima from a single interval LDL^T factorization."""
    size=len(b); R=b[0].parent()
    L=[[R(int(i==j)) for j in range(size)] for i in range(size)]
    D=[]; y=[]; values=[]
    distance=R(1)
    for i in range(size):
        for j in range(i):
            L[i][j]=(G[i][j]-sum((L[i][k]*D[k]*L[j][k]
                                      for k in range(j)),R(0)))/D[j]
        pivot=G[i][i]-sum((L[i][k]**2*D[k] for k in range(i)),R(0))
        assert pivot.is_finite() and pivot>0, ('nonpositive pivot',i+1)
        D.append(pivot)
        y.append(b[i]-sum((L[i][k]*y[k] for k in range(i)),R(0)))
        distance=distance-y[i]**2/pivot
        assert distance.is_finite() and distance>0 and distance<1
        values.append(distance)
    coefficients=[]
    for n in range(1,size+1):
        c=[y[i]/D[i] for i in range(n)]
        for i in range(n-1,-1,-1):
            c[i]-=sum((L[j][i]*c[j] for j in range(i+1,n)),R(0))
        coefficients.append(c)
    return values, coefficients, D, L, y


@lru_cache(None)
def computation(limit, digits):
    R,G,b=gram_system(limit,numberdb.bits(digits,losing=WORKING_GUARD))
    values,coefficients,pivots,L,y=ldl_distances(G,b)
    for v in values:
        assert v.accuracy() >= numberdb.bits(digits)+8
    return R,G,b,values,coefficients,pivots,L,y


class NymanBeurlingDistances(numberdb.Generator):
    table='T641'
    parameters=('N',)
    type='R'
    digits=100
    rigour='proven'
    files=('generate.py','checks.py','argument-selection.md','derivation.md',
           'numerical-evidence.json','skill-evidence.json')

    def enumerate(self, limit=MAX_N):
        for N in range(1,limit+1):
            yield {'N':str(N)}

    def value(self, params, digits):
        N=int(params['N'])
        if N<1: raise ValueError('N must be positive')
        return computation(max(MAX_N,N),digits)[3][N-1]


if __name__=='__main__':
    if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
        os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
    generator=NymanBeurlingDistances()
    if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
        print(generator.publish(message='certify squared Nyman-Beurling distances for N=1,...,30'))
    else:
        report=generator.verify(sample=None)
        print(report)
        sys.exit(0 if report.ok else 1)