back to table · edit · history · where entries came from · files · download
6566 bytes, as of the version from 2026-10-26 13:28 (current). Recorded here, not run.
"""Real zeros of the Scorer function Gi -- numberdb.org/T516.
Reproduce with SageMath and numberdb (mpmath is included in SageMath):
sage -pip install numberdb
sage -python generate.py # compare all twelve stored entries
sage -python generate.py --check # local certification, no network
sage -python generate.py --fill # fill this draft, NUMBERDB_API_KEY set
This does not publish or offer a draft for review. Every computation establishes
both the enclosures and the indexing anew. mpmath only proposes narrow brackets;
Arb certifies them. The first twelve negative roots are enumerated from zero.
At 700 bits, a 900-term Maclaurin polynomial encloses Gi and its derivative on
[-16,0]. Each of the three residue-class tails is bounded geometrically, using
its first omitted term and its decreasing ratio. For derivative order j=0,1,2,
the ratios are respectively R^3/((k+3)(k+2)), R^3/((k+2)k), R^3/(k(k-1)).
For indexing, on [m-h,m+h], Gronwall applied to (y,y') with the maximum norm
gives M=(|y(m)|+|y'(m)|+h/pi)*exp(16h) as a bound for both |y| and |y'|.
Thus y(m)+[-hM,hM] and y'(m)+[-h(16M+1/pi),h(16M+1/pi)] enclose y and y'.
Each segment either excludes zero in y, or has derivative of fixed sign and
its endpoint signs determine whether it has zero or one root. Inconclusive
segments are bisected. No homogeneous Sturm theorem is assumed.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
import mpmath as mp
BITS = 700
TERMS = 900
RADIUS = 16
ROOT_RADIUS = QQ(1) / 10**120
RF = RealBallField(BITS)
def rational_decimal(s):
"""Parse a finite decimal exactly, without a binary floating conversion."""
s = str(s).lower()
mantissa, _, exponent = s.partition('e')
exponent = int(exponent or 0)
whole, dot, frac = mantissa.partition('.')
return QQ(int(whole + frac)) * QQ(10)**(exponent-len(frac))
@lru_cache(None)
def coefficients():
c = [RF(1)/(RF(3)**(QQ(7)/6)*RF(QQ(2)/3).gamma()),
RF(1)/(RF(3)**(QQ(5)/6)*RF(QQ(1)/3).gamma()),
-RF(1)/(2*RF.pi())]
for k in range(TERMS):
c.append(c[k]/((k+3)*(k+2)))
return c
@lru_cache(None)
def polynomial_and_tail(order):
c = coefficients()
def factor(k):
f = 1
for j in range(order):
f *= k-j
return f
p = [c[k]*factor(k) for k in range(order, TERMS)]
tail = RF(0)
for k in range(TERMS, TERMS+3):
if order == 0:
denominator = (k+3)*(k+2)
elif order == 1:
denominator = (k+2)*k
elif order == 2:
denominator = k*(k-1)
else:
raise ValueError('supported derivative orders: 0, 1, 2')
q = RF(RADIUS)**3/denominator
assert q < 1
tail += abs(c[k])*factor(k)*RF(RADIUS)**(k-order)/(1-q)
assert tail.is_finite() and tail < RF(10)**-160
return p, tail
def gi(x, order=0):
x = RF(x)
assert abs(x) <= RADIUS
p, tail = polynomial_and_tail(order)
result = RF(0)
for c in reversed(p):
result = result*x+c
result = result.add_error(tail.upper())
assert result.is_finite()
return result
def sign(y):
if y > 0:
return 1
if y < 0:
return -1
return 0
def local_enclosures(a, b):
m, h = (a+b)/2, (b-a)/2
y, yp = gi(m), gi(m, 1)
M = (abs(y)+abs(yp)+RF(h)/RF.pi())*(RF(RADIUS)*h).exp()
ybox = y.add_error((RF(h)*M).upper())
dbox = yp.add_error((RF(h)*(RADIUS*M+1/RF.pi())).upper())
assert ybox.is_finite() and dbox.is_finite()
return ybox, dbox
@lru_cache(None)
def certify():
"""Return twelve rational brackets and a full finite-interval root count."""
stats = {'excluded_by_value': 0, 'excluded_by_monotonicity': 0,
'root_segments': 0, 'subdivisions': 0}
coarse = []
def visit(a, b, depth=0):
assert depth < 30
ybox, dbox = local_enclosures(a,b)
if sign(ybox):
stats['excluded_by_value'] += 1
return
if sign(dbox):
sa, sb = sign(gi(a)), sign(gi(b))
assert sa and sb
if sa != sb:
coarse.append((a,b))
stats['root_segments'] += 1
else:
stats['excluded_by_monotonicity'] += 1
return
stats['subdivisions'] += 1
m = (a+b)/2
visit(m,b,depth+1)
visit(a,m,depth+1)
for i in range(1, 16*16+1):
a,b = -QQ(i)/16,-QQ(i-1)/16
visit(a,b)
if len(coarse) >= 12:
assert len(coarse) == 12
stats['left_endpoint'] = str(a)
break
assert len(coarse) == 12
brackets=[]
for a,b in coarse:
with mp.workdps(170):
ma,mb = mp.mpf(str(a.numerator()))/int(a.denominator()),mp.mpf(str(b.numerator()))/int(b.denominator())
candidate=mp.findroot(mp.scorergi,(ma,mb))
mid=rational_decimal(mp.nstr(candidate,160))
lo,hi=mid-ROOT_RADIUS,mid+ROOT_RADIUS
assert a < lo < hi < b
assert sign(gi(lo))*sign(gi(hi)) == -1
assert sign(local_enclosures(lo,hi)[1])
brackets.append((lo,hi))
assert all(brackets[i][0] > brackets[i+1][1] for i in range(11))
return brackets,stats
class ScorerGiZeros(numberdb.Generator):
files = ('generate.py', 'checks.py')
table = 'T516'
parameters = ('n',)
type = 'R'
digits = 100
rigour = 'proven'
def enumerate(self):
for n in range(1,13):
yield {'n':n}
def value(self,params,digits):
if digits > 100:
raise ValueError('This certificate is configured for 100 digits.')
n=int(params['n'])
if not 1 <= n <= 12:
raise ValueError('The certified prefix has twelve roots.')
a,b=certify()[0][n-1]
return RF((a+b)/2).add_error(RF((b-a)/2).upper())
if __name__ == '__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
g=ScorerGiZeros()
if '--check' in sys.argv:
print(certify()[1])
for p in g.enumerate():
print(p,g.value(p,100))
elif '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(g.publish(message='Certify the first twelve signed real zeros, including their indices'))
else:
report=g.verify(sample=None)
print(report)
sys.exit(0 if report.ok else 1)