back to table · edit · history · where entries came from · files · download
1966 bytes, as of the version from 2026-10-28 17:16 (current). Recorded here, not run.
"""Certify Gram arguments for numberdb.org/T586.
Root routine reused from the live T583 generator, fetched 2026-10-28.
Analytic log-gamma fixes the phase continuously. Newton only proposes a
candidate; strict endpoint signs certify its enclosure. See certification.md.
"""
import os
import sys
import re
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
WORKING_GUARD = 160
ROOT_GUARD_DIGITS = 10
def phase(t):
"""Continuous Riemann-Siegel phase, using analytic log-gamma."""
R = t.parent()
C = ComplexBallField(R.precision())
return C(QQ(1)/4, t/2).log_gamma().imag() - t*R.pi().log()/2
def phase_derivative(t):
R = t.parent()
C = ComplexBallField(R.precision())
return (C(QQ(1)/4, t/2).psi().real() - R.pi().log())/2
def gram_point(n, digits=100):
if not isinstance(n, int) or n < -1:
raise ValueError('The phase index must be an integer at least -1')
R = RealBallField(numberdb.bits(digits, losing=WORKING_GUARD))
assert phase_derivative(R(7)) > 0
assert phase(R(7)) < -R.pi()
target = n*R.pi()
t = R(10 + 6*(n+1))
# No claim that a Newton candidate is an enclosure of the root.
for _ in range(20):
t_next = R((t - (phase(t)-target)/phase_derivative(t)).mid())
if t_next == t:
break
t = t_next
center = t.mid().exact_rational()
radius = QQ(1)/10**(digits+ROOT_GUARD_DIGITS)
left, right = center-radius, center+radius
assert left > 7
low_value = phase(R(left)) - target
high_value = phase(R(right)) - target
assert low_value.is_finite() and high_value.is_finite()
assert low_value < 0 < high_value, ('failed endpoint certificate', n)
ball = R(left).union(R(right))
assert ball.is_finite() and phase_derivative(ball) > 0
assert ball.accuracy() >= numberdb.bits(digits)
return ball