back to table · edit · history · where entries came from · files · download
6868 bytes, as of the version from 2026-10-25 08:29 (current). Recorded here, not run.
"""Tammes spherical-code separations -- numberdb.org/T479.
Run with SageMath and the attached sources.json:
sage -pip install numberdb
sage -python generate.py
sage -python checks.py
NUMBERDB_API_KEY=... sage -python generate.py --fill
The default verifies every stored entry. --fill writes entries and attachments
on the existing draft; it does not publish it or offer it for review.
The source snapshot was fetched on 2026-10-25. It uses Cohn's 3-dimensional
codes for 7..100 and Hardin, Sloane and Smith for 4..6. Published minimal
polynomials give 100-digit algebraic-cosine rows. Other rows retain ten
significant digits, with the printed source cosine enclosed by +/- 1e-12.
The 768-bit field gives at least 100 reliable decimal digits after root
isolation and conversion; checks.py measures the margin. Rounded coordinates
are independent checks, never the basis for extra digits or optimality claims.
"""
import os
import sys
import json
import re
from math import comb
from functools import lru_cache
from pathlib import Path
import numberdb.sage as numberdb
from sage.rings.rational_field import QQ
from sage.rings.real_arb import RealBallField
from numberdb._write import to_text
DATA = json.loads(Path(__file__).with_name('sources.json').read_text())
WORKING_BITS = 768
ROOT_STEPS = 400
def rational(s):
from decimal import Decimal
a,b=Decimal(str(s)).as_integer_ratio()
return QQ(a)/QQ(b)
def coefficients(text):
out={}
for term in re.findall(r'[+-]?[^+-]+',text):
if 'x' not in term: c,k=int(term),0
else:
left,right=term.split('x')
c=int(left.rstrip('*')) if left.rstrip('*') not in ('','+','-') else (-1 if left=='-' else 1)
k=int(right[1:]) if right else 1
out[k]=c
return [QQ(out.get(i,0)) for i in range(max(out)+1)]
def evaluate(cs,x):
ans=0
for c in reversed(cs): ans=ans*x+c
return ans
def isolate(cs,center):
"""Exact bisection, preceded by a rational uniqueness certificate."""
radius=QQ(1)/10**10
lo,hi=center-radius,center+radius
# Shift p(center+y) exactly to control its derivative without cancellation.
shifted=[sum(cs[j]*comb(j,k)*center**(j-k) for j in range(k,len(cs))) for k in range(len(cs))]
tail=sum(abs(shifted[k])*k*radius**(k-1) for k in range(2,len(cs)))
assert abs(shifted[1])>tail, 'initial interval does not certify uniqueness'
flo,fhi=evaluate(cs,lo),evaluate(cs,hi)
assert flo*fhi<0, 'source interval does not bracket a root'
for _ in range(ROOT_STEPS):
mid=(lo+hi)/2
fm=evaluate(cs,mid)
if fm==0:return mid,mid
if flo*fm<0: hi=mid
else: lo,flo=mid,fm
assert evaluate(cs,lo)*evaluate(cs,hi)<0
return lo,hi
@lru_cache(None)
def cosine(n):
R=RealBallField(WORKING_BITS)
if n==4:return R(-QQ(1)/3)
if n in (5,6):return R(0)
if n==9:return R(QQ(1)/3)
row=DATA['rows'][str(n)]
c=rational(row['cosine'])
if 'polynomial' in row:
a,b=isolate(coefficients(row['polynomial']),c)
return R((a+b)/2,(b-a)/2)
return R(c,QQ(1)/10**12)
def entry_comment(n):
if n==4:return r'The regular tetrahedron has $\cos(\pi\theta_4/180)=-1/3$ and $d_4=\sqrt{8/3}$. Optimality is proven CITE{FejesToth}.'
if n in (5,6):return r'$\theta_n=90$ and $d_n=\sqrt2$. Optimality is proven CITE{Schutte} CITE{FejesToth}.'
row=DATA['rows'][str(n)]
parts=[]
if 'polynomial' in row:
poly=row['polynomial'].replace('*','')
# At a 1365-pixel desktop viewport, the 50-character cutoff keeps
# every inline polynomial within the numbers pane. Longer algebraic
# descriptions remain in the attached source data.
if len(poly)<=50:
poly=re.sub(r'\^(\d+)',r'^{\1}',poly)
c=row['cosine']
parts.append(r'With $c_n=\cos(\pi\theta_n/180)$, $c_n$ is the unique zero of $'+poly+r'$ in $('+c+r'-10^{-10},'+c+r'+10^{-10})$ CITE{Cohn}.')
else:parts.append(r'The cosine $c_n=\cos(\pi\theta_n/180)$ is algebraic. Its polynomial and isolating interval are specified in the attached sources.json and generate.py CITE{Cohn}.')
if n in (7,8,9):parts.append('Optimality is proven CITE{Schutte}.')
elif n in (10,11):parts.append('Optimality is proven CITE{Danzer}.')
elif n==12:parts.append('Optimality is proven CITE{FejesToth}.')
elif n==13:parts.append('Optimality is proven CITE{Musin13}.')
elif n==14:parts.append('Optimality is proven CITE{Musin14}.')
elif n==24:parts.append('Optimality is proven CITE{Robinson}.')
elif n==15:parts.append('Cohn marks this configuration as optimal, citing the 2026 preprint CITE{Tammes15}; that proof claim is distinguished here from the established published cases.')
return ' '.join(parts)
class Tammes(numberdb.Generator):
table='T479'
parameters=('n','normalisation')
type='R'
digits=100
rigour='heuristic (agreement-checked)'
files=('generate.py','checks.py','check.py','sources.json','checks-result.json','source-notes.md')
def enumerate(self):
for n in range(4,101):
for normalisation in ('angle','chord'):
yield {'n':n,'normalisation':normalisation}
def value(self,params,digits):
n=int(params['n']);norm=params['normalisation']
R=RealBallField(WORKING_BITS);c=cosine(n)
if norm=='angle':
value=90 if n in (5,6) else c.arccos()*180/R.pi()
elif norm=='chord':value=(2-2*c).sqrt()
else:raise ValueError(norm)
precise=n<=6 or 'polynomial' in DATA['rows'][str(n)]
if not precise:
# A ball straddling a rounding boundary can make to_text shorten
# several digits. Round its midpoint, then certify the WHOLE input
# ball fits the decimal's stated +/- one-last-place interval.
from decimal import Decimal, localcontext
mid=QQ(value.mid())
with localcontext() as context:
context.prec=80
text=format(Decimal(int(mid.numerator()))/Decimal(int(mid.denominator())),'.10g')
ulp=QQ(10)**Decimal(text).as_tuple().exponent
assert abs(value-R(rational(text)))<R(ulp)
value=text
out={'number':value,'digits':100 if precise else 10}
note=entry_comment(n)
if note:out['comment']=note
return out
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
os.environ['NUMBERDB_API_KEY']=sys.stdin.read().strip()
g=Tammes()
if '--fill' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1':
print(g.publish(message='Fill checked Tammes angle and chord snapshot; retain draft status'))
else:
report=g.verify(sample=None)
print(report)
raise SystemExit(0 if report.ok else 1)