back to table · edit · history · where entries came from · files · download
4344 bytes, as of the version from 2026-10-30 14:29. Recorded here, not run.
"""Eccentricities at axisymmetric bifurcations of Maclaurin spheroids -- numberdb.org/T627.
Run with SageMath (or passagemath) and the numberdb package:
sage -pip install numberdb # once
sage -python generate.py # verify all stored entries
sage -python generate.py --fill-draft # fill the draft; requires API key
The --fill-draft option sends entries, never publishes or offers the table.
The default enumeration is k=1,...,10. See argument-selection.md and
certification.md for selection and proofs. No rounded source value selects
or initializes a root. Exact rational polynomials certify global uniqueness;
Arb signs bisect (0,1), then enclose e=1/sqrt(1+xi^2).
At 100 digits the working field is 500 bits; bisection retains at least 120
relative decimal digits in every e. The fixed 160-bit guard is also checked
against a 700-bit run. Failure of any certificate stops the calculation.
"""
import os
import sys
from functools import lru_cache
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.real_arb import RealBallField
WORKING_GUARD = 160
MIN_BITS = 500
@lru_cache(None)
def polynomials(k):
k = int(k)
if k < 1:
raise ValueError('The branch index must be positive.')
ring = PolynomialRing(QQ, 'x')
x = ring.gen()
p = [ring(1), x]
r = [ring(0), ring(1)]
ell = 2*k+2
for n in range(1, ell):
p.append(((2*n+1)*x*p[n] + n*p[n-1])/QQ(n+1))
r.append(((2*n+1)*x*r[n] + n*r[n-1])/QQ(n+1))
p, r = p[ell], r[ell]
# d(q/p)/dx = -1/((1+x^2)p^2), q=p atan(1/x)-r.
assert (1+x*x)*(r*p.derivative()-r.derivative()*p)-p*p == -1
assert all(c >= 0 for c in p.list()) and p[0] > 0
a, b = p*p+x*x, p*r+x
numerator = (1+x*x)*(b*a.derivative()-b.derivative()*a)-a*a
# H'=numerator/((1+x^2)A^2). Its numerator has one positive zero.
assert numerator[0] < 0
assert all(c >= 0 for c in numerator.list()[1:])
assert numerator.leading_coefficient() > 0
assert b.degree() < a.degree()
return p, r, a, b, numerator
def h(k, x):
_, _, a, b, _ = polynomials(k)
if x == 0:
return x.parent().pi()/2
return (1/x).arctan() - b(x)/a(x)
def enclose(k, digits=100, working_bits=None):
field = RealBallField(working_bits or max(MIN_BITS, numberdb.bits(digits, losing=WORKING_GUARD)))
lo, hi = QQ(0), QQ(1)
assert h(k, field(lo)) > 0 and h(k, field(hi)) < 0
# H(0)>0; H decreases then increases to zero. The root is unique on x>0.
steps = numberdb.bits(digits+22)
for _ in range(steps):
mid = (lo+hi)/2
v = h(k, field(mid))
if v > 0:
lo = mid
elif v < 0:
hi = mid
else:
raise ArithmeticError('Increase working precision: unresolved sign.')
xi = field((lo+hi)/2).add_error(field((hi-lo)/2))
e = 1/(1+xi*xi).sqrt()
assert e.is_finite() and e > 0 and e < 1
assert field(e.rad()) < field(10)**(-digits-18)
return e, lo, hi
class MaclaurinEccentricities(numberdb.Generator):
table = 'T627'
parameters = ('k',)
type = 'R'
digits = 100
rigour = 'proven'
files = ('generate.py', 'checks.py', 'argument-selection.md', 'certification.md',
'numerical-evidence.json', 'skill-evidence.json')
def enumerate(self, limit=10):
for k in range(1, int(limit)+1):
yield {'k': k}
def value(self, params, digits):
return enclose(int(params['k']), digits)[0]
def client_from_stdin():
if os.environ.get('NUMBERDB_KEY_FROM_STDIN') == '1':
key = sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isupper():
key = key.split('=',1)[1].strip().strip('\'"')
return numberdb.Client(api_key=key)
return None
if __name__ == '__main__':
client = client_from_stdin()
generator = MaclaurinEccentricities()
if '--fill-draft' in sys.argv or os.environ.get('NUMBERDB_PUBLISH') == '1':
print(generator.publish(client=client, message='Certify the first ten axisymmetric Maclaurin eccentricities at 100 digits.'))
else:
report = generator.verify(sample=None, client=client)
print(report)
sys.exit(0 if report.ok else 1)