Hilbert 12: closure attempt

Download original file

This is the original source text. Reading it does not run or verify the program.

1 #!/usr/bin/env python3
2 """Pure cubic H12 closing audit and exact-rational interval enclosure.
3 
4 Author: Komorebi P.
5 Uses only Python's standard library for the interval computation.
6 It DOES NOT prove the gamma value algebraic or its equality to U_alg.
7 Integer endpoints, explicit Taylor remainders, and an infinite-series tail
8 bound distinguish this calculation from floating-point recognition.
9 """
10 from __future__ import annotations
11 from dataclasses import dataclass
12 from fractions import Fraction
13 from math import factorial, isqrt
14 from collections import Counter
15 from pathlib import Path
16 import argparse
17 import json
18 import time
19 
20 BITS = 224
21 SCALE = 1 << BITS
22 
23 
24 def ceildiv(a: int, b: int) -> int:
25     if b <= 0:
26         raise ValueError('positive divisor required')
27     return -((-a) // b)
28 
29 
30 @dataclass(frozen=True)
31 class RI:
32     """Closed real interval [lo/SCALE, hi/SCALE]."""
33     lo: int
34     hi: int
35 
36     def __post_init__(self):
37         if self.lo > self.hi:
38             raise ValueError('reversed interval')
39 
40     @staticmethod
41     def of(x) -> 'RI':
42         if isinstance(x, RI):
43             return x
44         if isinstance(x, float):
45             raise TypeError('floats are not permitted')
46         q = Fraction(x)
47         return RI((q.numerator*SCALE)//q.denominator,
48                   ceildiv(q.numerator*SCALE, q.denominator))
49 
50     def __add__(self, x):
51         x=RI.of(x)
52         return RI(self.lo+x.lo, self.hi+x.hi)
53     __radd__=__add__
54 
55     def __neg__(self):
56         return RI(-self.hi, -self.lo)
57 
58     def __sub__(self,x):
59         return self + (-RI.of(x))
60 
61     def __rsub__(self,x):
62         return RI.of(x) + (-self)
63 
64     def __mul__(self,x):
65         x=RI.of(x)
66         vals=[self.lo*x.lo,self.lo*x.hi,self.hi*x.lo,self.hi*x.hi]
67         return RI(min(vals)//SCALE,ceildiv(max(vals),SCALE))
68     __rmul__=__mul__
69 
70     def reciprocal(self):
71         if self.lo <= 0 <= self.hi:
72             raise ZeroDivisionError('interval contains zero')
73         if self.hi < 0:
74             return -(-self).reciprocal()
75         return RI((SCALE*SCALE)//self.hi, ceildiv(SCALE*SCALE,self.lo))
76 
77     def __truediv__(self,x):
78         return self*RI.of(x).reciprocal()
79 
80     def square(self):
81         high=max(self.lo*self.lo,self.hi*self.hi)
82         low=0 if self.lo<=0<=self.hi else min(self.lo*self.lo,self.hi*self.hi)
83         return RI(low//SCALE,ceildiv(high,SCALE))
84 
85     def sqrt(self):
86         if self.lo < 0:
87             raise ValueError('negative lower bound')
88         low=isqrt(self.lo*SCALE); high=isqrt(self.hi*SCALE)
89         if high*high < self.hi*SCALE:
90             high+=1
91         return RI(low,high)
92 
93     def maxabs(self):
94         return max(abs(self.lo),abs(self.hi))
95 
96     def contains(self,x):
97         x=Fraction(x)
98         return self.lo*x.denominator <= x.numerator*SCALE <= self.hi*x.denominator
99 
100     def width(self):
101         return Fraction(self.hi-self.lo,SCALE)
102 
103 
104 @dataclass(frozen=True)
105 class CI:
106     re: RI
107     im: RI
108 
109     @staticmethod
110     def of(x) -> 'CI':
111         if isinstance(x, CI):
112             return x
113         return CI(RI.of(x), RI.of(0))
114 
115     def __add__(self,x):
116         x=CI.of(x)
117         return CI(self.re+x.re,self.im+x.im)
118     __radd__=__add__
119 
120     def __neg__(self):
121         return CI(-self.re,-self.im)
122 
123     def __sub__(self,x):
124         return self+(-CI.of(x))
125 
126     def __rsub__(self,x):
127         return CI.of(x)+(-self)
128 
129     def __mul__(self,x):
130         x=CI.of(x)
131         return CI(self.re*x.re-self.im*x.im,self.re*x.im+self.im*x.re)
132     __rmul__=__mul__
133 
134     def conjugate(self):
135         return CI(self.re,-self.im)
136 
137     def norm2(self):
138         return self.re.square()+self.im.square()
139 
140     def reciprocal(self):
141         d=self.norm2()
142         if d.lo <= 0:
143             raise ZeroDivisionError('complex box not separated from zero')
144         return CI(self.re/d,(-self.im)/d)
145 
146     def __truediv__(self,x):
147         return self*CI.of(x).reciprocal()
148 
149     def __rtruediv__(self,x):
150         return CI.of(x)*self.reciprocal()
151 
152     def __pow__(self,n:int):
153         if n<0:
154             return self.reciprocal()**(-n)
155         result=CI.of(1);x=self
156         while n:
157             if n&1:
158                 result=result*x
159             x=x*x;n//=2
160         return result
161 
162     def exp(self,terms:int=72):
163         # Reduce the complex modulus using the l1 upper bound to <= 1/2.
164         z=self; squarings=0
165         while z.re.maxabs()+z.im.maxabs() > SCALE//2:
166             z=z/2;squarings+=1
167         term=CI.of(1);total=term
168         for n in range(1,terms+1):
169             term=term*z/n
170             total=total+term
171         # |tail| <= exp(1/2)*(1/2)^(N+1)/(N+1)! < 2/(N+1)!
172         # The looser final bound is exact and sufficient at BITS=224.
173         error=ceildiv(2*SCALE,factorial(terms+1))
174         pad=RI(-error,error)
175         total=CI(total.re+pad,total.im+pad)
176         for _ in range(squarings):
177             total=total*total
178         return total
179 
180 
181 def atan_reciprocal(q:int,terms:int=130):
182     """Rational enclosure from the alternating arctan series."""
183     total=Fraction(0)
184     for k in range(terms):
185         total += Fraction((-1)**k,(2*k+1)*q**(2*k+1))
186     next_term=Fraction((-1)**terms,(2*terms+1)*q**(2*terms+1))
187     return min(total,total+next_term),max(total,total+next_term)
188 
189 
190 def pi_interval():
191     a,b=atan_reciprocal(5);c,d=atan_reciprocal(239)
192     low=16*a-4*d; high=16*b-4*c
193     return RI(RI.of(low).lo,RI.of(high).hi)
194 
195 
196 def icbrt(n:int):
197     lo=0;hi=1 << ((n.bit_length()+2)//3)
198     if hi**3 == n:
199         return hi
200     while hi-lo>1:
201         mid=(lo+hi)//2
202         if mid**3 <= n:lo=mid
203         else:hi=mid
204     return lo
205 
206 
207 def data():
208     sq3=RI.of(3).sqrt()
209     omega=CI(RI.of(Fraction(-1,2)),sq3/2)
210     ai=icbrt(7*SCALE**3)
211     alpha=RI(ai,ai+1)
212     beta=CI.of(alpha)*omega**2
213     pi=pi_interval()
214     twopii=CI(RI.of(0),2*pi)
215     return omega,beta,twopii
216 
217 
218 def decimal_out(i:RI,places:int=45):
219     """Outward decimal rendering, also without binary floating point."""
220     scale=10**places
221     low=(i.lo*scale)//SCALE
222     high=ceildiv(i.hi*scale,SCALE)
223     def fmt(n):
224         sign='-' if n<0 else ''
225         a=abs(n)
226         return f'{sign}{a//scale}.{a%scale:0{places}d}'
227     return '['+fmt(low)+', '+fmt(high)+']'
228 
229 
230 def box_out(z:CI,places=45):
231     return {'real':decimal_out(z.re,places),'imag':decimal_out(z.im,places)}
232 
233 
234 def algebraic_candidate(omega,beta,twopii):
235     # This is the retained, independently specified Gauss branch.
236     zeta=(twopii*Fraction(3,7)).exp()
237     chi=[None,CI.of(1),omega,omega**2,omega**2,omega,CI.of(1)]
238     gauss=CI.of(0);power=CI.of(1)
239     for a in range(1,7):
240         power=power*zeta
241         gauss=gauss+chi[a]*power
242     xi=gauss/beta
243     nu=(xi+beta/xi-beta**2)/3
244     v=-beta**2*nu-2*nu**2
245     Y=(4*beta**2+beta-14)*v**2+(11*beta-23)*v-beta**2+3*beta-3
246     Z=((-2*beta**2+7*beta-8)*v**2+(-10*beta**2+11*beta+17)*v)/3+7-beta-beta**2
247     u=(Y-(1+2*omega)*Z)/2
248     return u,xi,v
249 
250 
251 def gamma_value(omega,beta,twopii,terms=3000,evaluation_point=Fraction(1,3)):
252     tau=(beta**2+2*beta)/15
253     t=-tau;s=-(beta+2)/15
254     p=(twopii*t).exp();q=(twopii*s).exp()
255     r=Fraction(471,500);small=Fraction(1,2)
256     # Exact comparisons of squared absolute values with rational envelopes.
257     assert p.norm2().hi*r.denominator**2 < r.numerator**2*SCALE
258     assert q.norm2().hi*small.denominator**2 < small.numerator**2*SCALE
259     pn=CI.of(1);qn=CI.of(1);summ=CI.of(0)
260     if evaluation_point == Fraction(1,3):
261         roots=[CI.of(1),omega,omega**2]
262     elif evaluation_point == Fraction(1,4):
263         ii=CI(RI.of(0),RI.of(1))
264         roots=[CI.of(1),ii,CI.of(-1),-ii]
265     else:
266         raise ValueError('This certificate supports only the two stated rational points')
267     period=len(roots)
268     for n in range(1,terms+1):
269         pn=pn*p;qn=qn*q
270         if n%5:
271             summ=summ+(pn*roots[n%period]-qn*roots[-n%period])/(n*(1-pn)*(1-qn))
272     log_u=-5*summ
273     # Rigorous absolute tail for real z=1/3, excluding multiples of five.
274     tail=10*r**(terms+1)/((terms+1)*(1-r)**2*(1-small))
275     err=ceildiv(tail.numerator*SCALE,tail.denominator)
276     pad=RI(-err,err)
277     log_u=CI(log_u.re+pad,log_u.im+pad)
278     return log_u.exp(),log_u,tail,p,q
279 
280 
281 def action(x):
282     k,l,j=x
283     return ((4*k+5*l+j)%7,(4*l+5*j)%7,4*j%7)
284 
285 
286 def transport_audit():
287     pts=[(k,l,j) for k in range(7) for l in range(7) for j in range(7)]
288     ix={p:n for n,p in enumerate(pts)}
289     perm=[ix[action(p)] for p in pts]
290     assert perm[0]==0 and sorted(perm)==list(range(343))
291     visited=set();orbits=[]
292     for start in range(343):
293         if start in visited:continue
294         cyc=[];j=start
295         while j not in visited:
296             visited.add(j);cyc.append(j);j=perm[j]
297         assert j==start;orbits.append(cyc)
298     counts=Counter(map(len,orbits))
299     assert counts=={1:1,3:2,21:16}
300     # Check the full 21-step weighted return on every standard basis position.
301     # The abstract f_m are exact positive rationals, deliberately non-arithmetic.
302     def f(m,j):return Fraction((m+2)*(j+3)+5,m+j+7)
303     for j in range(343):
304         current=j; coefficient=Fraction(1)
305         for m in range(21,0,-1):
306             coefficient*=f(m,current)/f(m-1,perm[current])
307             current=perm[current]
308         assert current==j
309         assert coefficient==f(21,j)/f(0,j)
310     # Exact action of the fixed-point projector for arbitrary rational weights.
311     for j in range(343):
312         out_at_zero = (perm[j]==0)
313         in_at_zero=(j==0)
314         assert out_at_zero==in_at_zero
315     # Fourier orthogonality at the fixed coordinate is exact modulo Phi_7.
316     for x in pts:
317         counts_res=[0]*7
318         for a in pts:
319             counts_res[sum(ai*xi for ai,xi in zip(a,x))%7]+=1
320         if x==(0,0,0):assert counts_res==[343,0,0,0,0,0,0]
321         else:assert counts_res==[49]*7
322     return {'orbit_lengths':dict(counts),'weighted_21_step_checks':343,
323             'fixed_projector_checks':343,'Fourier_orthogonality_checks':343}
324 
325 
326 def interval_selftest():
327     tests=0
328     samples=[Fraction(-7,3),Fraction(-1,100),Fraction(0),Fraction(1,7),Fraction(5,2)]
329     for a in samples:
330         for b in samples:
331             ia=RI.of(a);ib=RI.of(b)
332             assert (ia+ib).contains(a+b)
333             assert (ia*ib).contains(a*b)
334             tests+=2
335             if b:
336                 assert (ia/ib).contains(a/b);tests+=1
337     # exp(0) must enclose 1 exactly.
338     e=CI.of(0).exp(); assert e.re.contains(1) and e.im.contains(0)
339     return tests+1
340 
341 
342 def main():
343     parser=argparse.ArgumentParser(description=__doc__)
344     parser.add_argument('--terms',type=int,default=3000)
345     parser.add_argument('--json',type=Path)
346     args=parser.parse_args()
347     if args.terms<1000:raise ValueError('Use at least 1000 terms for this closing audit')
348     start=time.monotonic()
349     result={'status':'analytic_arithmetic_equality_NOT_proved','interval_bits':BITS,
350             'selftests':interval_selftest(),'transport':transport_audit()}
351     omega,beta,twopii=data()
352     ua,xi,v=algebraic_candidate(omega,beta,twopii)
353     ug,logu,tail,p,q=gamma_value(omega,beta,twopii,args.terms)
354     ratio=ug/ua;cube=ratio**3;diff=ratio-1
355     bound=diff.re.maxabs()**2+diff.im.maxabs()**2
356     digits=0
357     for k in range(1,100):
358         if bound*10**(2*k)<SCALE*SCALE:digits=k
359         else:break
360     assert digits>=40, f'Insufficient ratio enclosure: only {digits} digits'
361     assert cube.re.lo > SCALE//2
362     assert diff.norm2().hi*4 < SCALE
363     # An enclosure containing zero does not certify exact equality.
364     assert diff.re.contains(0) and diff.im.contains(0)
365     quarter,_,_,_,_=gamma_value(omega,beta,twopii,args.terms,Fraction(1,4))
366     quarter_distance=(quarter/ua-1).norm2()
367     assert quarter_distance.lo*10000 > SCALE
368     result.update({'negative_control_quarter_distance':'> 1/100',
369                    'negative_control_quarter_value':box_out(quarter),
370                    'terms':args.terms,'U_gamma':box_out(ug),'U_alg':box_out(ua),
371                    'ratio':box_out(ratio),'ratio_cubed':box_out(cube),
372                    'certified_ratio_distance_upper':f'10^-{digits}',
373                    'certified_ratio_cubed_real_lower':'> 1/2',
374                    'certified_branch_disk':'|U_gamma/U_alg-1| < 1/2',
375                    'tail_exact_factored':f'10*(471/500)^({args.terms+1})/(({args.terms+1})*(29/500)^2*(1/2))',
376                    'tail_dyadic_upper_numerator':ceildiv(tail.numerator*SCALE,tail.denominator),
377                    'ratio_integer_endpoints':{'real':[ratio.re.lo,ratio.re.hi],
378                                               'imag':[ratio.im.lo,ratio.im.hi],
379                                               'common_denominator':SCALE},
380                    'zero_imaginary_part_proved':False,
381                    'elapsed_seconds':round(time.monotonic()-start,3)})
382     # Do not print the several-thousand-digit rational tail in the summary.
383     report={k:v for k,v in result.items() if not k.startswith('tail_')}
384     print(json.dumps(report,indent=2))
385     if args.json:
386         args.json.write_text(json.dumps(result,indent=2)+'\n')
387 
388 if __name__=='__main__':
389     main()

Back to source controls

Back to Hilbert 12: closure attempt