back to table · edit · history · where entries came from · files · download
10082 bytes, as of the version from 2026-10-31 10:42 (current). Recorded here, not run.
"""Independent certificates for T636; sage -python checks.py [--stored].
Download generate.py, T635.json and T79.json beside this file. By default,
check every generated decimal enclosure. --stored reads the draft (key in
NUMBERDB_API_KEY); NUMBERDB_KEY_FROM_STDIN=1 accepts the key from stdin.
Root certificates use rational sign brackets or a Taylor/Rouche inequality,
not QQbar's root algorithm. Snapshot comparisons, Laurent cancellation,
external polynomial examples and known-answer controls are also included.
"""
import os
import sys
import json
import re
from pathlib import Path
from fractions import Fraction
import numberdb.sage as numberdb
from numberdb._write import to_text
from sage.rings.integer_ring import ZZ
from sage.rings.rational_field import QQ
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.complex_arb import ComplexBallField
from sage.rings.real_mpfi import RealIntervalField
from sage.rings.complex_interval_field import ComplexIntervalField
from sage.arith.misc import binomial
from generate import CLASSES, MAX_DEGREE, FaberZeros, faber, coefficients
ROOT=Path(__file__).resolve().parent
QX=PolynomialRing(QQ,'X'); X=QX.gen()
F=ComplexBallField(768)
def scalar_bounds(s):
s=s.strip()
if s.startswith('['):
lo,hi=s[1:-1].split(','); return rational(lo),rational(hi)
q=rational(s)
if '.' not in s and 'e' not in s.lower(): return q,q
mantissa,*exponent=s.lower().split('e')
decimals=len(mantissa.split('.')[1]) if '.' in mantissa else 0
unit=QQ(10)**((int(exponent[0]) if exponent else 0)-decimals)
return q-unit,q+unit
def rational(s):
f=Fraction(s.strip()); return QQ(f.numerator)/QQ(f.denominator)
def rectangle(s):
if ' + i * ' in s:
a,b=s.split(' + i * ')
else: a,b=s,'0'
return scalar_bounds(a)+scalar_bounds(b)
def finite(z):
return z.real().is_finite() and z.imag().is_finite()
def sign_bracket(p,lo,hi):
return lo<hi and p(lo)*p(hi)<0
def certificate(p,box):
"""Prove a root in this rectangle without using algebraic root isolation."""
lo,hi,il,ih=box
p=p//p.gcd(p.derivative())
if il==ih==0:
if lo==hi:
assert p(lo)==0; return 'exact'
assert sign_bracket(p,lo,hi)
return 'real-sign'
if lo==hi:
# P(a+i*y) has rational real and imaginary polynomial parts.
shifted=p(X+lo)
real=QX([shifted[j]*(-1)**(j//2) if j%2==0 else 0 for j in range(p.degree()+1)])
imag=QX([shifted[j]*(-1)**((j-1)//2) if j%2 else 0 for j in range(p.degree()+1)])
common=real.gcd(imag); common=common//common.gcd(common.derivative())
assert common.degree()>0 and sign_bracket(common,il,ih)
return 'vertical-sign'
assert lo<hi and il<ih
c=F((lo+hi)/2,(il+ih)/2)
# Refine the disk center independently. A rectangle's midpoint can be
# farther from its root than the radius of its inscribed disk.
for _ in range(2):
c -= p(c)/p.derivative()(c)
assert finite(c)
RIF=RealIntervalField(768)
real=QQ(RIF(c.real()).center()); imag=QQ(RIF(c.imag()).center())
radius=min(real-lo,hi-real,imag-il,ih-imag)/2
assert radius>0, 'Newton center left the claimed rectangle'
c=F(real,imag)
# |P'(c)|r > |P(c)| + sum_{j>=2}|P^(j)(c)/j!|r^j.
terms=[]; derivative=p; factorial=ZZ(1)
for j in range(p.degree()+1):
if j:
derivative=derivative.derivative(); factorial*=j
term=derivative(c)/factorial
assert finite(term)
terms.append(abs(term))
linear=terms[1]*radius
remainder=terms[0]+sum((terms[j]*radius**j for j in range(2,len(terms))),F(0).real())
assert linear>remainder, ('Rouche inequality failed',box)
return 'rouche'
def rejected(function):
try: function()
except (AssertionError,ValueError): return
raise AssertionError('negative control was accepted')
def controls():
assert certificate(X**2-2,rectangle('1.41421356237309504880'))=='real-sign'
assert certificate(X**2-4,rectangle('2'))=='exact'
assert certificate(X**2+2,rectangle('0 + i * 1.41421356237309504880'))=='vertical-sign'
assert certificate(X**2-2*X+2,rectangle('1.0000000000 + i * 1.0000000000'))=='rouche'
rejected(lambda:certificate(X**2-2,rectangle('1.41421356237309504890')))
rejected(lambda:certificate(X**2-2*X+2,rectangle('1.0000000000 + i * 1.0100000000')))
assert not finite(F('nan'))
print('Known-answer and deliberately wrong-enclosure controls passed.',flush=True)
def principal_part(p,g,n):
# Sparse Laurent multiplication, independent of the Faber recurrence.
c=coefficients(g,n-1); t={-1:QQ(1)}
t.update({j:QQ(c[j]) for j in range(1,n)})
powers=[{0:QQ(1)}]
for d in range(1,n+1):
a={}
for i,x in powers[-1].items():
for j,y in t.items():
if i+j<=n-d:
a[i+j]=a.get(i+j,QQ(0))+x*y
powers.append(a)
for k in range(-n,1):
assert sum((p[d]*powers[d].get(k,QQ(0)) for d in range(n+1)),QQ(0)) == (1 if k==-n else 0)
def flatten(doc):
block=doc['Numbers']
if isinstance(block,list): return block
out=[]
for g,ns in block.items():
for n,ks in ns.items():
for k,value in ks.items():
row=dict(value) if isinstance(value,dict) else {'number':str(value)}
row['params']={'g':g,'n':n,'k':k}; out.append(row)
return out
def run(stored=False):
controls()
companion=json.loads((ROOT/'T635.json').read_text())['Numbers']
jdata=json.loads((ROOT/'T79.json').read_text())['Numbers']
assert all(coefficients('1A',9)[m]==ZZ(jdata[str(m)]) for m in range(1,10))
assert list(coefficients('13A',9)) == [0,12,28,66,132,258,468,843,1428,2406]
print('Nine positive j coefficients agree with T79; nine 13A coefficients agree with OEIS A034319.',flush=True)
# External examples from BKM pp.24 and 26 use j and T_2A+104.
for g,shift,examples in [('1A',744,[X**2-1488*X+159768,
X**3-2232*X**2+1069956*X-36866976]),
('2A',104,[X**2-208*X+2072,
X**3-312*X**2+19332*X-49568])]:
for n,example in enumerate(examples,2):
assert faber(g,n)==example(X+shift)
print('Four external polynomial examples agree after the stated shifts.',flush=True)
gen=FaberZeros()
if stored:
rows=flatten(numberdb.table('T636'))
else:
rows=[]
for params in gen.enumerate():
row=gen.value(params,100)
rows.append(dict(row,params=params,number=to_text(row['number'],100)))
groups={}; counts={}; minimum=10000; maxlen=0
for row in rows:
params=row['params']; g=params['g']; n=int(params['n']); k=int(params['k'])
groups.setdefault((g,n),[]).append((k,row))
assert list(groups)==[(g,n) for g in CLASSES for n in range(2,MAX_DEGREE+1)]
certificates={}; rational_rows=[]; nonreal=0
for (g,n),items in groups.items():
p=QX(companion[g][str(n)])
assert p==faber(g,n)
principal_part(p,g,n)
assert p.gcd(p.derivative()).degree()==0, 'update checker for repeated roots'
assert len(items)==n
boxes=[]
for k,row in items:
assert k==len(boxes)+1
box=rectangle(str(row['number'])); boxes.append(box)
kind=certificate(p,box)
certificates[kind]=certificates.get(kind,0)+1
assert row['comment']=='Multiplicity 1.'
if kind=='exact': rational_rows.append([g,n,k,str(row['number'])])
nonreal+=bool(box[2]!=0 or box[3]!=0)
value=gen.value({'g':g,'n':n,'k':k},100)['number']
if hasattr(value,'real'):
parts=[value.real(),value.imag()]
for part in parts:
if hasattr(part,'accuracy') and part!=0:
minimum=min(minimum,int(part.accuracy()))
# The full generated ball must lie inside the written interval.
RIF=RealIntervalField(768); CIF=ComplexIntervalField(768)
interval=CIF(value)
for part,(lo,hi) in zip([interval.real(),interval.imag()],[(box[0],box[1]),(box[2],box[3])]):
assert QQ(part.lower())>=lo and QQ(part.upper())<=hi
maxlen=max(maxlen,len(str(row['number'])))
# Disjoint rectangles plus n existence certificates account for all n roots.
for i,a in enumerate(boxes):
for b in boxes[i+1:]:
assert a[1]<b[0] or b[1]<a[0] or a[3]<b[2] or b[3]<a[2]
# Independent ordering: strict real separation, or exactly conjugate boxes.
for a,b in zip(boxes,boxes[1:]):
if a[1]<b[0]: continue
assert a[:2]==b[:2] and a[2]==-b[3] and a[3]==-b[2] and a[3]<b[2]
# Vieta and conjugation, evaluated on the actual written intervals.
vals=[CIF(RealIntervalField(768)(a,b),RealIntervalField(768)(c,d)) for a,b,c,d in boxes]
assert 0 in sum(vals,CIF(0)).real() and 0 in sum(vals,CIF(0)).imag()
prod=CIF(1)
for z in vals: prod*=z
assert 0 in (prod-((-1)**n)*p[0]).real() and 0 in prod.imag()
for a in boxes:
assert any(a[:2]==b[:2] and a[2]==-b[3] and a[3]==-b[2] for b in boxes)
counts[g]=counts.get(g,0)+len(items)
print('Certified',g,n,'roots',len(items),flush=True)
print('RESULT_JSON',json.dumps({'entries':len(rows),'polynomials':len(groups),'counts':counts,
'certificates':certificates,'nonreal':nonreal,'rational_rows':rational_rows,
'minimum_component_bits':minimum,'longest_written_value':maxlen,'stored':stored}))
if not stored: print('ENTRIES_JSON',json.dumps(rows))
if __name__=='__main__':
if os.environ.get('NUMBERDB_KEY_FROM_STDIN')=='1':
key=sys.stdin.read().strip()
if '=' in key and key.split('=',1)[0].isidentifier(): key=key.split('=',1)[1].strip().strip('\"\'')
numberdb.configure(api_key=key)
run(stored='--stored' in sys.argv or os.environ.get('NUMBERDB_PUBLISH')=='1')