back to table · edit · history · where entries came from · files · download
4937 bytes, as of the version from 2026-11-09 23:44 (current). Recorded here, not run.
"""Log-sine-cosine integrals -- numberdb.org/T699.
Run with SageMath (or passagemath with numberdb):
sage -pip install numberdb
sage -python generate.py # verify the stored draft/table
sage -python generate.py --submit # write entries with your API key
The generator uses Arb polynomial arithmetic and analytic logarithmic
endpoint moments, with a positive-coefficient bound for the entire omitted
tail. See certification.md. All selection logic is in this file. Read
argument-selection.md for the finite coverage rationale and skill-manifest.json
for the archived instructions used in this build. No cached numerical data
are required. Submission fills the existing draft; it never offers or publishes.
"""
import os
import sys
from functools import lru_cache
from math import comb, factorial
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
WORKING_GUARD_BITS = 192
TAIL_GUARD_DEGREE = 32
@lru_cache(maxsize=None)
def series_data(t_string, digits):
"""Power series in w=(v/a)^2, where a=pi*t <= pi/2."""
t = QQ(t_string)
assert 0 < t <= QQ(1)/2
bits = numberdb.bits(digits, losing=WORKING_GUARD_BITS)
degree = bits + TAIL_GUARD_DEGREE
R = RealBallField(bits)
P = PolynomialRing(R, 'w')
a = R.pi()*R(t)
ratio = R(t/2)**2
c = [R(0)] + [R(2*j).zeta()*ratio**j/j for j in range(1, degree+1)]
sine = P([a.log()] + [-x for x in c[1:]])
cosine = P([R(2).log()] + [-(4**j-1)*c[j] for j in range(1, degree+1)])
return R, P, a, degree, sine, cosine
@lru_cache(maxsize=None)
def powers(t_string, digits, which, exponent):
R, P, a, degree, sine, cosine = series_data(t_string, digits)
if exponent == 0:
return P.one()
base = sine if which == 'sine' else cosine
return (powers(t_string, digits, which, exponent-1)*base).truncate(degree+1)
@lru_cache(maxsize=None)
def segment(m, n, t_string, digits):
"""Negative integral from 0 to pi*t, t<=1/2, with a proved tail bound."""
p, q = m-1, n-1
R, P, a, degree, sine, cosine = series_data(t_string, digits)
value = R(0)
cosine_power = powers(t_string, digits, 'cosine', q)
for i in range(p+1):
poly = (powers(t_string, digits, 'sine', p-i)*cosine_power).truncate(degree+1)
# integral_0^1 u^(2*j) log(u)^i du = (-1)^i i!/(2*j+1)^(i+1).
moment = sum((coefficient/R(2*j+1)**(i+1)
for j, coefficient in enumerate(poly.list())), R(0))
value += comb(p, i)*(-1)**i*factorial(i)*moment
value *= -a
# Majorants evaluated at w=2, valid uniformly for 0<t<=1/2.
A = -2*(R(7)/8).log()
B = 2*R(2).log()
majorant = sum((R(comb(p,i)*factorial(i))*(abs(a.log())+A)**(p-i)
for i in range(p+1)), R(0))*(R(2).log()+B)**q
error = a*R(2)**(-degree-1)*majorant
value = value.add_error(error.upper())
assert value.is_finite() and value.accuracy() >= numberdb.bits(digits+20)
return value
@lru_cache(maxsize=None)
def compute(m, n, t_string, digits=100):
t = QQ(t_string)
assert m >= 1 and n >= 1 and 0 < t <= 1
if t <= QQ(1)/2:
return segment(m,n,str(t),digits)
full = segment(m,n,'1/2',digits) + segment(n,m,'1/2',digits)
if t == 1:
return full
return full - segment(n,m,str(1-t),digits)
class LogSineCosine(numberdb.Generator):
table = 'T699'
parameters = ('m','n','t')
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py','argument-selection.md','certification.md','checks.py','skill-manifest.json')
def enumerate(self):
rows = []
for total in range(4,11):
for m in range(2,total-1):
n = total-m
rows.append((m,n,QQ(1)/2))
if m <= n:
rows.append((m,n,QQ(1)))
rows.extend([(3,2,QQ(1)/6),(2,3,QQ(1)/6),
(3,2,QQ(1)/3),(3,2,QQ(5)/6)])
for m,n,t in sorted(rows,key=lambda row:(row[0]+row[1],row[0],row[1],row[2])):
yield {'m':str(m),'n':str(n),'t':str(t)}
def value(self,params,digits):
return compute(int(params['m']),int(params['n']),params['t'],digits)
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key = sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isidentifier():
key = key.split('=',1)[1].strip().strip('\"\'')
os.environ['NUMBERDB_API_KEY'] = key
generator = LogSineCosine()
if '--submit' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(message='Compute 48 certified mixed logarithmic moments with reproducible checks'))
else:
report = generator.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)