#!/usr/bin/env python3 """Exact rational certificates for the revised Borsuk manuscript. Run with Python 3. No third-party packages or input files are needed. Checks: main base <1.166374, small-parameter base <1.1673, quartic base <1.166474, and the stated method-specific floor sqrt(843)/25. Analytic lemmas are proved in the manuscript; this program verifies their displayed numerical certificates. """ from fractions import Fraction as F def verify_main_bound(): from fractions import Fraction as F from math import factorial def rat(s):return F(str(s)) def trim(p): while len(p)>1 and p[-1]==0:p.pop() return p def evaluate(p,x): v=F(0) for c in reversed(p):v=v*x+c return v def derivative(p):return [i*p[i] for i in range(1,len(p))] def remainder(p,q): p=trim(p.copy()) while len(p)>=len(q) and any(p): k=len(p)-len(q);r=p[-1]/q[-1] for i,c in enumerate(q):p[i+k]-=r*c trim(p) return p def sign(x):return (x>0)-(x<0) def changes(seq): signs=[sign(y) for y in seq if y] return sum(a!=b for a,b in zip(signs,signs[1:])) def sturm_counts(poly): ss=[poly,derivative(poly)] while len(ss[-1])>1: rem=[-v for v in remainder(ss[-2],ss[-1])] if not any(rem):break scale=abs(rem[-1]);ss.append([v/scale for v in rem]) return changes([p[0] for p in ss]),changes([p[-1] for p in ss]) gamma,eta,lam,m0,c0,c1,c2,theta6,theta10,qb,qcut,qr,qt,cv,A0,xi=map(rat,['0.46943966493022', '0.03147922506626', '0.83422581962954', '0.84031353613959', '-0.03861305150538', '0.29647789431113', '-0.96425074771268', '0.01246697666460', '0.00826314934702', '0.00363340143362', '3.75487290159014', '0.77810686755460', '0.16218889900094', '0.88617587397029', '0.11733062741006', '0.00322745373409']) d=16*cv-6-4*lam alpha=6*d/(6*A0+d)**2;beta0=36*A0*A0/(6*A0+d)**2 fcoef4=(16*cv*cv-4*lam*lam)/3;fcoef2=2+2*lam*lam-8*cv*cv k4=2*gamma-gamma*gamma*(1-eta)*fcoef4-xi;k3=-2*gamma;k6=-F(8,17)*gamma*gamma*eta co=[F(0)]*21 co[0]=-c0;co[1]=-c1;co[2]=-c2;co[3]=k3;co[4]=k4;co[6]=k6 co[10]-=theta6/945;co[12]+=theta6/10395 co[18]-=theta10/34459425;co[20]+=theta10/654729075 Qpoly=[F(0)]*8;Qpoly[3]=qb*qcut**4;Qpoly[5]=-2*qb*qcut*qcut;Qpoly[7]=qb fq=[F(0)]*15 for i,z in enumerate(Qpoly):fq[i+1]+=2*z;fq[i]-=2*z Qpoly[1]-=qr for i,zz in enumerate(Qpoly): for j,ww in enumerate(Qpoly):fq[i+j]-=zz*ww for i,z in enumerate(fq):co[i]+=qt*z st=sturm_counts(co) assert st[0]==st[1] and co[0]>0 and co[-1]>0 print('Degree-20 polynomial Sturm variations:',st) co_margin=co.copy();co_margin[0]-=F(1,1000000) st_margin=sturm_counts(co_margin) assert st_margin[0]==st_margin[1] and co_margin[0]>0 and co_margin[-1]>0 print('Polynomial exceeds 1e-6 on the nonnegative half-line:',st_margin) assert 00 and c1>0 and theta6>0 and theta10>0 and A0>0 assert d>0 and fcoef2<0 and alpha>0 and beta0>0 and fcoef4>0 f0=c0-2*gamma*gamma*(1-eta)*lam*lam-(1-eta)*beta0 f1=c1+c2+4*gamma*gamma*(1-eta)*lam f2=-gamma*gamma*(1-eta)*fcoef2;f3=F(18,17)*gamma*gamma*eta P0=f1+2*f2*m0+3*f3*m0*m0 Q0=f2*m0*m0+2*f3*m0**3-f0 kA=qt+eta+(1-eta)*alpha beta,rho,s,sigma,zeta,M,omega=map(rat,['0.7866529541', '0.8456140263', '0.1255274413', '0.6908074308', '0.6774479657', '1.3417650561', '0.0128557668']) q=sigma*sigma;u=s*s/q assert 0=0 and M>0 and omega>0 from math import comb # Coefficients in y=1-q. K_a[y^k]=a/(a-k)*(y^k-y^a). def flow_polys(n): p=[F((-1)**k*comb(n,k)) for k in range(n+1)];arr=[p] for r in range(2,9): ar=r*(2*r-1);pnew=[F(0)]*(ar+1) for k,z in enumerate(p): z*=F(ar,ar-k);pnew[k]+=z;pnew[ar]-=z p=pnew;arr.append(p) return arr entropy={} for n in (1,3,5): ys=[evaluate(p,1-q) for p in flow_polys(n)] assert all(y>0 for y in ys) entropy[n]=sum(u**r*ys[r-1]/r for r in range(1,9))+ys[-1]*u**9/(9*(1-u)) assert sigma**4<=F(3,5) def correction_polys(qpoly): p=[sum(qpoly[l]*comb(l,k)*(-1)**k for l in range(k,len(qpoly))) for k in range(len(qpoly))] arr=[p] for r in range(3,9): ar=r*(2*r-1);pn=[F(0)]*(ar+1) for k,z in enumerate(p): z*=F(ar,ar-k);pn[k]+=z;pn[ar]-=z p=pn;arr.append(p) return arr def correction_sum(qpoly): ys=[evaluate(p,1-q) for p in correction_polys(qpoly)] assert all(y>0 for y in ys) return sum(u**r*ys[r-2]/r for r in range(2,9))+ys[-1]*u**9/(9*(1-u)) U=correction_sum([F(0),F(0),F(1,2),F(-1),F(3,4),F(-3,10),F(1,20)]) V=correction_sum([F(0),F(0),F(0),F(0),F(1,4),F(-3,10),F(1,20)]) commgap=omega*xi-2*U assert commgap>0 A=(1+zeta*zeta)/(2*zeta);B=(1+zeta*zeta)/4 M0=1+A*(1-beta*rho)**2+B*beta*beta def bj(j):return s*s*q**(2*j-2)+A*(s*q**(j-1)-beta*rho**(2*j))**2-B*beta*beta def entropyj(j):return 2*j*(entropy[2*j-1] if j<3 else entropy[5]*q**(2*j-6)) def dj(j):return bj(j)/M+entropyj(j) gap1=2*omega*P0+24*(U*m0+V)-dj(1) gap2=-4*omega*kA-dj(2) gap3=-6*omega*qt-dj(3) gap4=-8*omega*qt-dj(4) gaps=[-dj(j) for j in range(5,51)] assert gap1>0 and gap2>0 and gap3>0 and gap4>0 and min(gaps)>0 assert q*q0 assert B*beta*beta/(2*M)>F(84,1000) x=M0/M-1+omega*Q0+6*U*m0*m0+12*V assert 00 for name,val in [('d',d),('alpha',alpha),('beta0',beta0),('kA',kA),('U',U),('V',V),('commutator absorption margin',commgap),('P0',P0),('Q0',Q0),('entropy n1',entropy[1]),('entropy n3',entropy[3]),('entropy n5',entropy[5]),('degree 1 margin',gap1),('degree 2 margin',gap2),('degree 3 margin',gap3),('degree 4 margin',gap4),('degrees 5..50 margin',min(gaps)),('tail margin',tailgap),('squared-base margin',basegap)]:print(name,float(val)) print('ALL RATIONAL CHECKS PASSED: C < 1.166374') """Small-parameter elementary certificate, using fractions only.""" from fractions import Fraction as F def verify_elementary(): m0=F(5,6) slope=F(27,16)*m0-F(3,8) intercept=F(27,32)*m0*m0-F(3,8)*m0-F(3,16)-slope*m0 assert slope==F(33,32) and intercept==-F(99,128) A0=F(1,8) ds=1-F(3,2)*A0 dc=F(3,4)*A0*A0 assert ds==F(13,16) and dc==F(3,256) assert ds*F(2,3)+dc>F(1,3) P=slope/ds;Q=(dc-intercept)/ds assert P==F(33,26) and Q==F(201,208) beta,rho,s,sigma=F(4,5),F(5,6),F(1,8),F(3,4) zeta,M,omega=F(2,3),F(4,3),F(1,120) A=(1+zeta*zeta)/(2*zeta); B=(1+zeta*zeta)/4 M0=1+A*(1-beta*rho)**2+B*beta*beta assert M0==F(3649,2700) u=s*s/(sigma*sigma); cbar=s*s/(1-u) assert cbar==F(9,560) def d(j): bj=s*s*sigma**(4*j-4)+A*(s*sigma**(2*j-2)-beta*rho**(2*j))**2-B*beta*beta return bj/M+2*j*cbar*sigma**(4*j-4) assert d(1)==F(306967,14515200) assert d(2)==-F(20591657753,300987187200) assert 2*omega*P-d(1)==F(1109,188697600)>0 assert -4*omega-d(2)==F(10558751513,300987187200)>0 assert sigma**80 and F2<0 and k4>0 c0=-k4*bb*bb;c1=-2*aa*bb*k4;c2=-k4*(aa*aa-2*bb) assert c1>0 assert alpha==F(594,625) and beta0==F(9,625) and k4==F(281,1200) assert aa==F(540,281) assert -2*aa*k4==-2*gamma and -c1==2*aa*bb*k4 and -c0==k4*bb*bb and -c2==k4*(aa*aa-2*bb) # k4*x^4-2*gamma*x^3-c0-c1*x-c2*x^2 = k4*(x^2-aa*x-bb)^2. f0=c0-2*gamma*gamma*lam*lam-beta0;f1=c1+c2+4*gamma*gamma*lam;f2=-gamma*gamma*F2 P=f1+2*f2*m0;Q=f2*m0*m0-f0 beta,rho,s,sigma,zeta,M,omega=map(rat,['0.7827996394', '0.8491703718', '0.1249486022', '0.7381553257', '0.6785218289', '1.3420283373', '0.0137606737']) q=sigma*sigma;u=s*s/q assert 0=0 and M>0 and omega>0 cplus=q*(sum(u**k/k for k in range(1,9))+u**9/(9*(1-u))) A=(1+zeta*zeta)/(2*zeta);B=(1+zeta*zeta)/4 M0=1+A*(1-beta*rho)**2+B*beta*beta bj=lambda j:s*s*q**(2*j-2)+A*(s*q**(j-1)-beta*rho**(2*j))**2-B*beta*beta dj=lambda j:bj(j)/M+2*j*cplus*q**(2*j-2) gap1=2*omega*P-dj(1);gap2=-4*omega*alpha-dj(2) gaps=[-dj(j) for j in range(3,51)] assert gap1>0 and gap2>0 and min(gaps)>0 assert q*q0 x=M0/M-1+omega*Q assert 00 for name,z in [('alpha',alpha),('beta0',beta0),('k4',k4),('aa',aa),('P',P),('Q',Q),('first margin',gap1),('second margin',gap2),('degrees3..50',min(gaps)),('infinite tail',tail),('squared-base margin',basegap)]:print(name,str(z) if name in ['alpha','beta0','k4','aa','P','Q'] else float(z)) print('ALL EXACT CHECKS PASSED: C<1.166474') def verify_method_floor(): from fractions import Fraction as F from itertools import product from math import comb DIM = 4 class Poly(dict): def __add__(self, other): if not isinstance(other, Poly): other = Poly({(0,) * DIM: F(other)}) ans = Poly(self) for key, value in other.items(): ans[key] = ans.get(key, F(0)) + value return Poly({key: value for key, value in ans.items() if value}) __radd__ = __add__ def __neg__(self): return Poly({key: -value for key, value in self.items()}) def __sub__(self, other): return self + (-other if isinstance(other, Poly) else -F(other)) def __rsub__(self, other): return -self + other def __mul__(self, other): if not isinstance(other, Poly): other = Poly({(0,) * DIM: F(other)}) ans = Poly() for key, value in self.items(): for key2, value2 in other.items(): key3 = tuple(x + y for x, y in zip(key, key2)) ans[key3] = ans.get(key3, F(0)) + value * value2 return Poly({key: value for key, value in ans.items() if value}) __rmul__ = __mul__ def __pow__(self, power): ans = Poly({(0,) * DIM: F(1)}) for _ in range(power): ans = ans * self return ans def indices(shape): return product(*(range(n) for n in shape)) def bernstein(poly): """Power coefficients to tensor-product Bernstein coefficients on [0,1]^4.""" shape = tuple(1 + max(key[a] for key in poly) for a in range(DIM)) ans = {key: poly.get(key, F(0)) for key in indices(shape)} for axis, size in enumerate(shape): old, ans = ans, {} for key in indices(shape): ans[key] = sum(old[key[:axis] + (i,) + key[axis + 1:]] * F(comb(key[axis], i), comb(size - 1, i)) for i in range(key[axis] + 1)) return shape, ans def split(data, axis): """Exact de Casteljau subdivision at the midpoint of coordinate axis.""" shape, coeffs = data left, right = {}, {} remaining = tuple(a for a in range(DIM) if a != axis) n = shape[axis] - 1 for rest in product(*(range(shape[a]) for a in remaining)): def key(i): out = [0] * DIM out[axis] = i for a, value in zip(remaining, rest): out[a] = value return tuple(out) row = [coeffs[key(i)] for i in range(n + 1)] left[key(0)], right[key(n)] = row[0], row[-1] for level in range(1, n + 1): row = [(row[i] + row[i + 1]) / 2 for i in range(len(row) - 1)] left[key(level)], right[key(n - level)] = row[0], row[-1] return (shape, left), (shape, right) def certify(target): r = Poly({(1,0,0,0):F(1)});s=Poly({(0,1,0,0):F(1)}) q = Poly({(0,0,1,0):F(1)});z=Poly({(0,0,0,1):F(1)}) h=F(7,11);coeff=[h/2,h/24,h/80] V=sum(coeff[j]*q**(2*j) for j in range(3)) C=sum(coeff[j]*q**j*r**(2*j+2) for j in range(3)) D=sum(coeff[j]*r**(4*j+4) for j in range(3)) W=sum((2*j+2)*coeff[j]*q**(2*j) for j in range(3)) a=1+z*z;d=2*(r*r+D)+h*z p=((2*z+a)*(1+s*s*V)*d-2*a*(r+s*C)**2)*(1+s*s*W)-2*target*z*d boundary=2*(sum(coeff[j]*(r**(2*j+2)-s*r*q**j)**2 for j in range(3)) +s*s*sum(coeff[j]*coeff[k]*(q**j*r**(2*k+2)-q**k*r**(2*j+2))**2 for j in range(3) for k in range(j+1,3)))*(1+s*s*W) remainder=p-boundary assert all(key[-1]>=1 for key in remainder) remainder=Poly({key[:-1]+(key[-1]-1,):value for key,value in remainder.items()}) # p = boundary + z * remainder, and boundary >= 0, z >= 0. stack = [(bernstein(p), bernstein(remainder), (0,) * DIM, 0)] visited = leaves = max_depth = 0 while stack: bp, br, depths, depth = stack.pop() visited += 1 if min(bp[1].values()) >= 0 or min(br[1].values()) >= 0: leaves += 1 max_depth = max(max_depth, depth) continue assert depth < 45, "Unexpected failure to certify positivity" axis = min(range(DIM), key=lambda a: depths[a]) lp, rp = split(bp, axis) lr, rr = split(br, axis) new_depths = depths[:axis] + (depths[axis] + 1,) + depths[axis + 1:] stack.append((lp, lr, new_depths, depth + 1)) stack.append((rp, rr, new_depths, depth + 1)) return visited, leaves, max_depth target=F(843,625) result=certify(target) assert result==(619,310,22) print("METHOD FLOOR CERTIFIED: squared base >=", target, "nodes/leaves/depth:", result) if __name__ == "__main__": verify_main_bound() result=verify_elementary() print("ELEMENTARY CERTIFICATE PASSED: base < 1.1673", result["base_squared"]) verify_quartic_bound() verify_method_floor() print("ALL FOUR EXACT NUMERICAL CERTIFICATES PASSED")