back to table · edit · history · where entries came from · files · download
4472 bytes, as of the version from 2026-10-28 17:16 (current). Recorded here, not run.
"""Certify Hardy Z stationary arguments for numberdb.org/T586.
The isolation and refinement routines are reused from the live attached
T585 generator, fetched 2026-10-28. This module computes roots afresh;
no stored decimal argument is used. See certification.md for the proof.
"""
import os
import sys
import json
import hashlib
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from sage.misc.cython import cython_import
from sage.rings.real_arb import RealBallField
from sage.rings.rational_field import QQ
COUNT = 100
COVER_BITS = 96
WORKING_GUARD = 160
@lru_cache(maxsize=1)
def bridge():
return cython_import(str(Path(__file__).with_name('hardy_arb.pyx')))
def sign(x):
if not x.is_finite():
return 0
return 1 if x > 0 else -1 if x < 0 else 0
def jet(a, b, prec, length=3):
R = RealBallField(prec)
x = R((QQ(a)+QQ(b))/2, (QQ(b)-QQ(a))/2)
result = bridge().jet(x, length, prec)
# The FLINT Hardy jet is real for real arguments. Nonfinite/wide balls
# are never used as successful certificates; they cause subdivision.
return [v.real() for v in result]
@lru_cache(maxsize=None)
def isolate(count=COUNT):
roots = []
leaves = []
origin_end = QQ(1)/1024
origin_second = 2*jet(0, origin_end, COVER_BITS)[2]
assert sign(origin_second) == 1
leaves.append(('0', str(origin_end), 'origin', 1))
def cover(a, b, depth=0):
coefficients = jet(a, b, COVER_BITS)
first_sign = sign(coefficients[1])
if first_sign:
leaves.append((str(a), str(b), 'exclude', first_sign))
return
second_sign = sign(coefficients[2])
if second_sign:
sa = sign(jet(a, a, COVER_BITS, 2)[1])
sb = sign(jet(b, b, COVER_BITS, 2)[1])
if sa and sb:
if sa != sb:
assert sb == second_sign and sa == -second_sign
roots.append((a,b))
leaves.append((str(a), str(b), 'root', second_sign))
else:
leaves.append((str(a), str(b), 'monotone-exclude', sa))
return
if depth >= 30:
raise ArithmeticError('Unresolved interval: %s, %s' % (a,b))
midpoint = (a+b)/2
cover(a,midpoint,depth+1)
cover(midpoint,b,depth+1)
a = origin_end
k = 0
while len(roots) < count:
k += 1
if k > 4000:
raise ArithmeticError('Increase the explicit height budget')
b = QQ(k)/4
cover(a,b)
a = b
assert len(roots) == count
# Check exact adjacency, not just the total length of the cover.
cursor = QQ(0)
totals = {}
for left,right,method,s in leaves:
assert QQ(left) == cursor and QQ(right) > QQ(left)
cursor = QQ(right)
totals[method] = totals.get(method,0)+1
assert cursor == a and totals['root'] == count
encoded = json.dumps(leaves,separators=(',',':')).encode()
summary = {'count':count,'endpoint':str(a),'cover_bits':COVER_BITS,
'leaf_counts':totals,'leaves_sha256':hashlib.sha256(encoded).hexdigest(),
'origin_second_derivative':str(origin_second),
'origin_endpoint':str(origin_end)}
return tuple(roots), summary, leaves
@lru_cache(maxsize=None)
def root(j, digits=100):
if not 1 <= j <= COUNT:
raise ValueError('This certified implementation covers j=1,...,100')
(a,b) = isolate()[0][j-1]
prec = numberdb.bits(digits,losing=WORKING_GUARD)
R = RealBallField(prec)
x = R((a+b)/2)
# All division here is in the real ball field. Midpoints are used only
# as candidates, never as a claim of exactness for the root.
for _ in range(18):
values = bridge().jet(x,3,prec)
f = values[1].real()
df = 2*values[2].real()
assert sign(df)
new = x-f/df
assert new.is_finite()
x = R(new.mid())
assert R(a) < x and x < R(b)
centre = QQ(x.mid())
radius = QQ(1) / (2**numberdb.bits(digits+15,losing=0))
left, right = centre-radius, centre+radius
assert a < left < right < b
fl = jet(left,left,prec,2)[1]
fr = jet(right,right,prec,2)[1]
second = 2*jet(left,right,prec)[2]
assert sign(fl)*sign(fr) == -1
assert sign(second) == sign(fr)
result = R(centre,radius)
assert result.is_finite() and result.accuracy() >= numberdb.bits(digits)+24
return result