back to table · edit · history · where entries came from · files · download
7947 bytes, as of the version from 2026-10-28 00:06 (current). Recorded here, not run.
"""Independent angular-integral and serialization checks for numberdb.org/T565.
Run: sage -python checks.py
For checks of the stored draft, set NUMBERDB_CHECK_STORED=1 and provide a key.
The Fourier determinant and Pfaffians use no Gamma factors and no generator
helpers. Their derivation and known controls are in independent-checks.md.
"""
import os
import sys
import json
from math import factorial, comb
from decimal import Decimal
from fractions import Fraction
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from numberdb._write import to_text
def binomial(n, k):
return QQ(comb(n, k)) if 0 <= k <= n else QQ(0)
def determinant(rows):
a = [list(map(QQ, row)) for row in rows]
result = QQ(1)
for i in range(len(a)):
pivot = next((j for j in range(i, len(a)) if a[j][i]), None)
if pivot is None:
return QQ(0)
if pivot != i:
a[i], a[pivot] = a[pivot], a[i]
result = -result
v = a[i][i]
result *= v
for r in range(i + 1, len(a)):
ratio = a[r][i] / v
for s in range(i + 1, len(a)):
a[r][s] -= ratio*a[i][s]
return result
def pfaffian(rows):
a = [list(map(QQ, row)) for row in rows]
assert len(a) % 2 == 0
result = QQ(1)
for i in range(0, len(a), 2):
pivot = next((j for j in range(i + 1, len(a)) if a[i][j]), None)
if pivot is None:
return QQ(0)
if pivot != i + 1:
a[i+1], a[pivot] = a[pivot], a[i+1]
for row in a:
row[i+1], row[pivot] = row[pivot], row[i+1]
result = -result
v = a[i][i+1]
result *= v
for r in range(i+2, len(a)):
for s in range(r+1, len(a)):
a[r][s] -= (a[i][r]*a[i+1][s] - a[i][s]*a[i+1][r])/v
a[s][r] = -a[r][s]
return result
def sign_power(integer):
return QQ(-1 if int(integer) % 2 else 1)
def skew_mode(a, b):
"""Integral sign(y-x) exp(i*a*x+i*b*y) dxdy, divided by i*pi.
a,b belong to the same integer or half-integer lattice.
"""
if a and b == -a:
return -QQ(4)/a
if a == 0 and b:
assert b.denominator() == 1
return -4*sign_power(b)/b
if b == 0 and a:
assert a.denominator() == 1
return 4*sign_power(a)/a
return QQ(0)
@lru_cache(None)
def angular_integral(beta, n, k):
"""Exact (r,p) from the defining integral, with value r*pi**p."""
if beta == 2:
a = [[binomial(2*k, k+i-j) for j in range(n)] for i in range(n)]
return QQ(factorial(n))*determinant(a), 0
if beta == 4:
a = [[(j-i)*binomial(4*k, 2*k+i+j-(2*n-1))
for j in range(2*n)] for i in range(2*n)]
return QQ(factorial(n))*pfaffian(a), 0
assert beta == 1
modes = [[(QQ(2*j - (n-1) + k - 2*t)/2, binomial(k,t))
for t in range(k+1)] for j in range(n)]
size = n + (n % 2)
a = [[QQ(0) for j in range(size)] for i in range(size)]
for i in range(n):
for j in range(i+1,n):
a[i][j] = sum((ca*cb*skew_mode(fa,fb)
for fa,ca in modes[i] for fb,cb in modes[j]), QQ(0))
a[j][i] = -a[i][j]
border_pi = 0
if n % 2:
border_pi = int(k % 2 == 0)
for i in range(n):
value = QQ(0)
for f,c in modes[i]:
if border_pi:
if f == 0:
value += 2*c
else:
assert f.denominator() == 2
sine = sign_power((int(2*f)-1)//2)
value += c*2*sine/f
a[i][n], a[n][i] = value, -value
m = n//2
phase = sign_power(m) if n % 2 else QQ(1)
return phase*QQ(factorial(n))/QQ(2**n)*pfaffian(a), m+border_pi-n
def controls():
assert determinant([[2,3],[5,7]]) == -1
assert pfaffian([[0,2,3,5],[-2,0,7,11],[-3,-7,0,13],[-5,-11,-13,0]]) == 28
assert pfaffian([[0,0,2,0],[0,0,0,3],[-2,0,0,0],[0,-3,0,0]]) == -6
assert angular_integral(1,1,0) == (QQ(1),0)
assert angular_integral(1,1,1) == (QQ(4),-1)
assert angular_integral(1,2,0) == (QQ(4),-1)
assert angular_integral(1,3,0) == (QQ(6),-1)
assert angular_integral(2,2,0) == (QQ(2),0)
assert angular_integral(4,2,0) == (QQ(6),0)
return 9
def decimal_interval(text):
decimal = Decimal(text)
center = Fraction(decimal)
unit = Fraction(10)**decimal.as_tuple().exponent
return QQ(center-unit), QQ(center+unit)
def contains_written(text, coefficient, power):
if '.' not in text and 'e' not in text.lower():
assert power == 0 and QQ(text) == coefficient
return
rb = RealBallField(640)
ball = rb(coefficient)*rb.pi()**power
assert ball.is_finite() and ball.accuracy() >= 620
low, high = decimal_interval(text)
assert low < QQ(ball.lower().exact_rational()) and QQ(ball.upper().exact_rational()) < high, text
def barnes_integer(m):
result = QQ(1)
for j in range(1,m-1):
result *= factorial(j)
return result
def main():
from generate import MorrisIntegrals, exact_form
total_controls = controls()
generator = MorrisIntegrals()
entries, precisions = [], []
for params in generator.enumerate():
beta,n,k = (int(params[a]) for a in generator.parameters)
independent = angular_integral(beta,n,k)
generated = exact_form(beta,n,k)
assert generated == independent, (params,generated,independent)
value = generator.value(params,100)
text = to_text(value,100)
contains_written(text,*independent)
if independent[1]:
assert value.is_finite()
precisions.append(int(value.accuracy()))
if beta == 2:
g=barnes_integer
assert independent[0]/QQ(factorial(n)) == g(n+1)*g(n+2*k+1)*g(k+1)**2/(g(n+k+1)**2*g(2*k+1))
if k == 1:
assert independent == (QQ(factorial(n+1)),0)
if k == 1:
assert independent == angular_integral(beta,n+1,0)
entries.append({'params':params,'number':text,'rational_coefficient':str(independent[0]),'pi_power':independent[1]})
if os.environ.get('NUMBERDB_CHECK_STORED') == '1':
document = numberdb.table('T565')
stored = document['Numbers']
count=0
for beta, ns in stored.items():
for n, ks in ns.items():
for k, text in ks.items():
if isinstance(text,dict): text=text['number']
q,p=angular_integral(int(beta),int(n),int(k))
contains_written(str(text),q,p)
count+=1
assert count == 150
report=generator.verify(sample=None)
assert report.ok
print('Stored entries:',count,'independent interval inclusions; generator verify:',report)
result={'table':'T565','entries':len(entries),'exact_rows':sum(e['pi_power']==0 for e in entries),
'irrational_rows':sum(e['pi_power']!=0 for e in entries),'independent_integral_checks':len(entries),
'known_controls':total_controls,'working_bits':numberdb.bits(100,losing=96),
'minimum_retained_bits':min(precisions),'independent_enclosure_bits':640,
'longest_serialized':max(len(e['number']) for e in entries),
'longest_at':max(entries,key=lambda e:len(e['number']))['params'],
'serialized_values_bytes':sum(len(e['number']) for e in entries)}
print('EVIDENCE_JSON='+json.dumps(result,sort_keys=True))
print('ENTRIES_JSON='+json.dumps(entries))
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key=sys.stdin.read().strip()
if '=' in key:key=key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY']=key
main()