back to table · edit · history · where entries came from · files · download
4399 bytes, as of the version from 2026-11-11 07:57 (current). Recorded here, not run.
"""T713: Generalized log-sinh integrals, 100 significant digits.
Usage: sage -pip install numberdb
sage -python generate.py # verify all stored entries
sage -python generate.py --publish # fill the authorized draft
Requires SageMath or passagemath with real/complex balls and polynomials.
All endpoint inputs are reconstructed exactly in RealBallField. The defining
improper integral is evaluated by integrating a finite polynomial in t^2 and
log(t), with an explicit uniform remainder; see computation.md. The series
comes from DLMF 4.36.1. No numerical quadrature is used by this generator.
Working precision is bits(digits+80); the polynomial cutoff is digits+50.
At 100 requested digits the whole 98-row selection is checked separately by
mpmath quadrature and low-order/alternating-binomial identities (checks.py).
The default enumeration is a noncanonical finite selection, documented in
argument-selection.md; the method supports the same endpoints at other orders.
Credentials and session provenance must be supplied by the caller.
"""
import math
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
ENDPOINTS = ('logphi', '2logphi', '2logsilver')
def endpoint(field, key):
a = ((field(1) + field(5).sqrt()) / 2).log()
if key == 'logphi':
return a
if key == '2logphi':
return 2 * a
if key == '2logsilver':
return 2 * (field(1) + field(2).sqrt()).log()
raise ValueError('Unknown exact endpoint')
@lru_cache(maxsize=None)
def prepare(key, digits):
R = RealBallField(numberdb.bits(digits + 80))
x = endpoint(R, key)
q = (x / (2 * R.pi()))**2
assert 0 < q < QQ(1)/10
N = digits + 50
coefficients = [x.log()]
qj = R(1)
for j in range(1, N + 1):
qj *= q
coefficients.append((-1)**(j+1) * R(2*j).zeta() * qj / j)
polynomial = PolynomialRing(R, 'y2')(coefficients)
# zeta(2j)<2 gives a uniform bound for the omitted B(x*y) terms.
epsilon = 2*q**(N+1) / ((N+1)*(1-q))
assert 2*q/(1-q) < 1
return R, x, polynomial, epsilon
@lru_cache(maxsize=None)
def power_coefficients(key, digits, power):
R, x, polynomial, epsilon = prepare(key, digits)
return (polynomial**power).list()
class LogSinh(numberdb.Generator):
table = 'T713'
parameters = ('n', 'k', 'x')
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
for n in range(2, 11):
for k in range(n - 1):
for key in ENDPOINTS:
if key == 'logphi' and not (k == 0 and n <= 9):
continue
yield {'n': n, 'k': k, 'x': key}
@lru_cache(maxsize=None)
def compute(self, n, k, key, digits):
n, k, digits = int(n), int(k), int(digits)
if n < 1 or k < 0 or k >= n:
raise ValueError('Require n>=1 and 0<=k<n')
p = n - 1 - k
R, x, polynomial, epsilon = prepare(key, digits)
total = R(0)
for h in range(p + 1):
coeffs = power_coefficients(key, digits, p-h)
moment = sum((c / R(k+2*j+1)**(h+1)
for j, c in enumerate(coeffs)), R(0))
total += math.comb(p, h)*(-1)**h*math.factorial(h)*moment
answer = -x**(k+1) * total
# |log(x*y)| <= |log(x)| - log(y), and |B|,|B_N|<1.
C = abs(x.log()) + 1
bound = sum((R(math.comb(p-1,h)*math.factorial(h)) *
C**(p-1-h) / R(k+1)**(h+1)
for h in range(p)), R(0)) if p else R(0)
error = p * epsilon * x**(k+1) * bound
answer = answer.add_error(error)
assert answer.is_finite()
assert answer.accuracy() > numberdb.bits(digits+10)
return answer
def value(self, params, digits):
return self.compute(int(params['n']), int(params['k']),
str(params['x']), int(digits))
if __name__ == '__main__':
import os
import sys
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
numberdb.configure(api_key=sys.stdin.read().strip())
g = LogSinh()
if '--publish' in sys.argv:
print(g.publish(removing=False))
else:
print(g.verify(sample=None))