back to table · edit · history · where entries came from · files · download
4011 bytes, as of the version from 2026-10-25 20:38 (current). Recorded here, not run.
"""Independent numerical check for Bernstein's constant, numberdb.org/T503.
Run: sage -python remez_check.py
Uses mpmath through Sage; it does not import the transcription generator.
Computes even minimax polynomials for |x| on [-1,1], solves the alternating
error equations, finds every critical point on [0,1], and extrapolates in
1/n^2. All arithmetic is fixed precision: this is numerical evidence only.
The exact degree-two solution p(x)=x^2+1/8, E_2=1/8 is the control.
"""
import numberdb.sage as numberdb
import mpmath as mp
import json
import time
def horner(a, x):
y=mp.mpf(0)
for c in reversed(a):
y=y*x+c
return y
def root(f, df, lo, hi, tol):
flo=f(lo); fhi=f(hi)
assert flo*fhi<0
x=(lo+hi)/2
for _ in range(500):
fx=f(x)
if abs(fx)<tol:
return x
if flo*fx>0:
lo=x;flo=fx
else:
hi=x;fhi=fx
if hi-lo<tol:
return (lo+hi)/2
dfx=df(x)
candidate=x-fx/dfx if dfx else (lo+hi)/2
x=candidate if lo<candidate<hi else (lo+hi)/2
raise ArithmeticError('root did not converge')
def minimax(n, dps):
with mp.workdps(dps):
tol=mp.mpf(10)**(-dps+60)
xs=[mp.sin(mp.pi*j/(2*(n+1))) for j in range(n+2)]
xs[0]=mp.mpf(0);xs[-1]=mp.mpf(1)
for iteration in range(30):
A=mp.matrix([[x**(2*k) for k in range(n+1)]+[mp.mpf((-1)**j)] for j,x in enumerate(xs)])
coeff=mp.lu_solve(A,mp.matrix(xs));a=list(coeff)[:-1];err=abs(coeff[-1])
da=[k*a[k] for k in range(1,n+1)]
dda=[k*da[k] for k in range(1,n)]
f=lambda x:x-horner(a,x*x)
df=lambda x:1-2*x*horner(da,x*x)
ddf=lambda x:-2*horner(da,x*x)-4*x*x*horner(dda,x*x)
zeros=[root(f,df,xs[j],xs[j+1],tol) for j in range(n+1)]
extrema=[root(df,ddf,zeros[j],zeros[j+1],tol) for j in range(n)]
newxs=[mp.mpf(0)]+extrema+[mp.mpf(1)]
spread=max(abs(abs(f(x))-err) for x in newxs)
if spread<mp.mpf(10)**(-dps+70):
# The polynomial df has n positive roots found in disjoint
# brackets and n sign changes in its coefficients, so it has
# no additional positive roots (Descartes' rule).
coeff_desc=[mp.mpf(1)]+[-2*(k+1)*a[k+1] for k in range(n)]
signs=[mp.sign(c) for c in coeff_desc if c]
assert sum(signs[j]!=signs[j+1] for j in range(len(signs)-1))==n
assert all(0<x<1 for x in extrema)
assert all(f(newxs[j])*f(newxs[j+1])<0 for j in range(n+1))
return 2*n*err,a,spread,iteration+1
xs=newxs
raise ArithmeticError('Remez did not converge')
def extrapolate(values,dps):
with mp.workdps(dps):
vals=[mp.mpf(x) for x in values]
xs=[mp.mpf(1)/(n*n) for n in range(1,len(vals)+1)]
for k in range(1,len(vals)):
vals=[(xs[i]*vals[i+1]-xs[i+k]*vals[i])/(xs[i]-xs[i+k]) for i in range(len(vals)-1)]
return vals[0]
def main():
dps=160
with mp.workdps(dps):
s,a,spread,it=minimax(1,dps)
assert abs(s-mp.mpf(1)/4)<mp.mpf('1e-65')
assert abs(a[0]-mp.mpf(1)/8)<mp.mpf('1e-65')
assert abs(a[1]-1)<mp.mpf('1e-65')
print('CONTROL degree 2: p=x^2+1/8, E=1/8: PASS',flush=True)
values=[];rows=[];started=time.time()
for n in range(1,53):
s,a,spread,it=minimax(n,dps)
values.append(mp.nstr(s,125))
row={'n':n,'2nE':mp.nstr(s,120),'spread':mp.nstr(spread,8),'iterations':it}
rows.append(row)
print(json.dumps(row),flush=True)
estimates={str(n):mp.nstr(extrapolate(values[:n],dps),100) for n in (8,16,24,32,40,48,50,52)}
print('RESULT '+json.dumps({'dps':dps,'mpmath_version':mp.__version__,'seconds':time.time()-started,'rows':rows,'extrapolations':estimates}),flush=True)
if __name__=='__main__':
main()