Exact program for the Euler worked example
Python 3.10 or later and its standard library. Exact integers and fractions.
Read the full derivation Download Python · Source and results bundle · Source identities
Run python -B euler_packets.py. Add --output my-results to save exact tables.
1 #!/usr/bin/env python3
2 """Exact controls for the extended Euler packet, polynomial and residual case.
3
4 Python 3.10+, standard library only. No network, floating arithmetic or model API.
5 By default this writes only a concise report to stdout. --output DIRECTORY saves
6 JSON and CSV files there, refusing to replace existing result files.
7 """
8 from __future__ import annotations
9
10 import argparse
11 import csv
12 from fractions import Fraction as F
13 import hashlib
14 from itertools import permutations
15 import json
16 from math import comb, factorial
17 from pathlib import Path
18 import platform
19
20
21 SECANT = (1, 1, 5, 61, 1385, 50521, 2702765, 199360981,
22 19391512145, 2404879675441, 370371188237525, 69348874393137901,
23 15514534163557086905, 4087072509293123892361,
24 1252259641403629865468285, 441543893249023104553682821)
25 QUARTIC = (1, 1, 69, 33661, 60376809, 288294050521, 3019098162602349,
26 60921822444067346581, 2159058013333667522020689,
27 125339574046311949415000577841,
28 11289082167259099068433198467575829,
29 1510335441937894173173702826484473600301)
30 REFERENCES = {
31 'secant': dict(id='A000364', offset=0,
32 url='https://oeis.org/A000364/b000364.txt', checked='2026-10-29'),
33 'quartic': dict(id='A211212', offset=0,
34 url='https://oeis.org/A211212', checked='2026-10-29'),
35 }
36
37
38 def choose(n, k):
39 return comb(n, k) if 0 <= k <= n else 0
40
41
42 def trim(a):
43 a = list(a)
44 while len(a) > 1 and not a[-1]:
45 a.pop()
46 return a or [F(0)]
47
48
49 def add(a, b, length=None):
50 length = max(len(a), len(b)) if length is None else length
51 return [F(a[i] if i < len(a) else 0) + F(b[i] if i < len(b) else 0)
52 for i in range(length)]
53
54
55 def scale(a, scalar):
56 return [F(scalar) * v for v in a]
57
58
59 def mul(a, b, length=None):
60 length = len(a)+len(b)-1 if length is None else length
61 result = [F(0)] * length
62 for i, x in enumerate(a[:length]):
63 if x:
64 for j, y in enumerate(b[:length-i]):
65 if y:
66 result[i+j] += x*y
67 return result
68
69
70 def shift(a, amount, length):
71 return ([F(0)]*amount + list(a))[:length]
72
73
74 def inverse(a, length):
75 if not a or not a[0]:
76 raise ValueError('A formal unit needs a nonzero constant coefficient')
77 b = [1/F(a[0])]
78 for n in range(1, length):
79 b.append(-sum(F(a[j])*b[n-j] for j in range(1, min(n+1, len(a))))/a[0])
80 return b
81
82
83 def power(a, exponent, length=None):
84 if exponent < 0:
85 if length is None:
86 raise ValueError('Negative powers require a precision')
87 return power(inverse(a, length), -exponent, length)
88 result = [F(1)]
89 for _ in range(exponent):
90 result = mul(result, a, length)
91 return result
92
93
94 def evaluate(a, value):
95 total = F(0)
96 for coefficient in reversed(a):
97 total = total*value+coefficient
98 return total
99
100
101 def translate(a, offset, scalar=1):
102 # The polynomial a(offset + scalar*x), by the ordinary binomial theorem.
103 result = [F(0)]*len(a)
104 for n, value in enumerate(a):
105 for j in range(n+1):
106 result[j] += value*comb(n,j)*offset**(n-j)*scalar**j
107 return trim(result)
108
109
110 def reciprocal_recurrence(spacing, parameter, count):
111 """P_N,n(r) from D*F' + r*D'*F = 0, in compressed EGF coordinates."""
112 if spacing < 1 or count < 1:
113 raise ValueError('Positive spacing and count are required')
114 r = F(parameter)
115 values = [F(1)]
116 for n in range(1, count):
117 values.append(sum((-1)**(k+1)*(
118 choose(spacing*n-1, spacing*k)
119 + r*choose(spacing*n-1, spacing*k-1))*values[n-k]
120 for k in range(1,n+1)))
121 return values
122
123
124 def reciprocal_binomial(spacing, parameter, count):
125 """Independent ordinary Cauchy products for (1+U)^(-r)."""
126 u = [F(0)] + [F((-1)**n, factorial(spacing*n)) for n in range(1,count)]
127 result = [F(1)]+[F(0)]*(count-1)
128 term = [F(1)]+[F(0)]*(count-1)
129 coefficient = F(1)
130 for k in range(1,count):
131 term = mul(term, u, count)
132 coefficient *= -F(parameter+k-1)/k
133 result = add(result, scale(term, coefficient), count)
134 return [result[n]*factorial(spacing*n) for n in range(count)]
135
136
137 def reciprocal_polynomials(spacing, count):
138 values = [[F(1)]]
139 for n in range(1,count):
140 row = [F(0)]*(n+1)
141 for k in range(1,n+1):
142 weight = [choose(spacing*n-1,spacing*k),
143 choose(spacing*n-1,spacing*k-1)]
144 row = add(row,scale(mul(weight,values[n-k]),(-1)**(k+1)),n+1)
145 values.append(trim(row))
146 return values
147
148
149 def alternating_count(spacing, size):
150 """Rank DP: a descent occurs unless the position is divisible by spacing."""
151 if size == 0:
152 return 1
153 row = [1]
154 for position in range(1,size):
155 if position % spacing == 0: # preceding value is smaller
156 row = [sum(row[:j]) for j in range(position+1)]
157 else:
158 row = [sum(row[j:]) for j in range(position+1)]
159 return sum(row)
160
161
162 def difference_coefficients(values):
163 row, result = list(values), []
164 while row:
165 result.append(row[0])
166 row = [y-x for x,y in zip(row,row[1:])]
167 return result
168
169
170 def transfer_row(spacing, offset, order):
171 unit = [F(comb(spacing,j)) for j in range(spacing+1)]
172 unit[0] -= 1
173 prefix = [F(comb(offset,j)) for j in range(offset+1)]
174 return trim(mul(prefix,power(unit,order)))
175
176
177 def as_integer(value):
178 value = F(value)
179 if value.denominator != 1:
180 raise AssertionError('Expected an integer, got '+str(value))
181 return value.numerator
182
183
184 def local_source(chart, degree):
185 """Paper IX source, indexed by ordinary even degree."""
186 if degree < 0 or degree % 2:
187 raise ValueError('An even nonnegative ordinary degree is required')
188 if chart == 2:
189 return F(3**(degree+2)-1,4*(degree+1)*(degree+2))
190 if chart < 1 or chart % 2 == 0:
191 raise ValueError('Use the dyadic chart 2 or a positive odd length')
192 return F(sum((-1)**a*(2*a+1)**degree for a in range(chart)))
193
194
195 def local_inverse_kernel(chart, count):
196 """EGF coefficients of the explicit inverse transfer unit."""
197 if chart == 2:
198 return [F(2**(4*r+1),(2*r+1)*(2*r+2)) for r in range(count)]
199 if chart < 1 or chart % 2 == 0:
200 raise ValueError('Use the dyadic chart 2 or a positive odd length')
201 return [F(1)]+[F(2**(2*r-1)*chart**(2*r)) for r in range(1,count)]
202
203
204 def local_kernel(chart, count):
205 """Invert the explicit unit by ordinary rational Cauchy products."""
206 ordinary=[v/factorial(2*r) for r,v in enumerate(local_inverse_kernel(chart,count))]
207 return [v*factorial(2*r) for r,v in enumerate(inverse(ordinary,count))]
208
209
210 def local_terms(chart, index):
211 """Terms reconstructing E_(2*index), including every factorial weight."""
212 return [comb(2*index,2*r)*coefficient*local_source(chart,2*(index-r))
213 for r,coefficient in enumerate(local_kernel(chart,index+1))]
214
215
216 def bernoulli_numbers(degree):
217 """Independent inversion of (exp(t)-1)/t, with B_1=-1/2."""
218 ordinary=inverse([F(1,factorial(j+1)) for j in range(degree+1)],degree+1)
219 return [value*factorial(j) for j,value in enumerate(ordinary)]
220
221
222 def valuation(value, prime):
223 """Exact rational p-valuation; None denotes the valuation of zero."""
224 value=F(value)
225 if value == 0:
226 return None
227 numerator,denominator=abs(value.numerator),value.denominator
228 result=0
229 while numerator % prime == 0:
230 numerator//=prime
231 result+=1
232 while denominator % prime == 0:
233 denominator//=prime
234 result-=1
235 return result
236
237
238 def rational_residue(value, modulus):
239 """Reduce only fractions whose reduced denominator is a unit."""
240 value=F(value)
241 return value.numerator*pow(value.denominator,-1,modulus) % modulus
242
243
244 def local_row_record(chart, index):
245 terms=local_terms(chart,index)
246 exact=sum(terms)
247 partial=F(0)
248 rows=[]
249 for depth,term in enumerate(terms):
250 partial+=term
251 error=exact-partial
252 rows.append(dict(depth=depth,term=str(term),partial=str(partial),
253 error=str(error),valuation=valuation(error,chart)))
254 return dict(chart=chart,ordinary_degree=2*index,exact=str(exact),rows=rows)
255
256
257 def composition_count(spacing, blocks):
258 """Inclusion-exclusion over the cuts between descending blocks."""
259 if blocks == 0:
260 return 1
261 def compositions(total):
262 if total == 0:
263 yield ()
264 else:
265 for first in range(1,total+1):
266 for rest in compositions(total-first):
267 yield (first,)+rest
268 result=0
269 for parts in compositions(blocks):
270 denominator=1
271 for part in parts:
272 denominator*=factorial(spacing*part)
273 result+=(-1)**(blocks-len(parts))*factorial(spacing*blocks)//denominator
274 return result
275
276
277 def toeplitz_matrix(coefficients, size):
278 """Leading block of multiplication in ordinary coefficient coordinates."""
279 return [[F(coefficients[n-j]) if j <= n else F(0)
280 for j in range(size)] for n in range(size)]
281
282
283 def coefficient_matrix(kernel, size):
284 """Leading block in even exponential coefficient coordinates."""
285 return [[F(comb(2*n,2*j))*kernel[n-j] if j <= n else F(0)
286 for j in range(size)] for n in range(size)]
287
288
289 def matrix_product(a, b):
290 if not a or not b or any(len(row) != len(b) for row in a):
291 raise ValueError('Compatible nonempty matrices are required')
292 if any(len(row) != len(b[0]) for row in b):
293 raise ValueError('The right matrix must be rectangular')
294 return [[sum((x*y for x,y in zip(row,column)),F(0))
295 for column in zip(*b)] for row in a]
296
297
298 def matrix_vector(a, vector):
299 if any(len(row) != len(vector) for row in a):
300 raise ValueError('Matrix and vector dimensions differ')
301 return [sum((x*y for x,y in zip(row,vector)),F(0)) for row in a]
302
303
304 def pascal_matrix(size, inverse=False):
305 return [[F(choose(n,j)*((-1)**(n-j) if inverse and j <= n else 1))
306 for j in range(size)] for n in range(size)]
307
308
309 def congruential_packet(spacing, address, count):
310 """Normalized positive-section Euler coefficients C_(N,j)(n)."""
311 if spacing < 1 or not 0 <= address < spacing:
312 raise ValueError('Require N>=1 and 0<=address<N')
313 denominator=[F(factorial(address),factorial(spacing*n+address)) for n in range(count)]
314 return [value*factorial(spacing*n) for n,value in enumerate(inverse(denominator,count))]
315
316
317 def fibonacci_polynomials(count):
318 """B_k=F_(k+1), ascending ordinary polynomial coefficients."""
319 basis=[]
320 previous,current=[F(0)],[F(1)]
321 for _ in range(count):
322 basis.append(current)
323 previous,current=current,trim(add([F(0)]+current,previous))
324 return basis
325
326
327 def monic_coordinates(polynomial, basis):
328 """Descending elimination, independent of the ballot coefficient formula."""
329 remainder=list(polynomial)+[F(0)]*(len(basis)-len(polynomial))
330 coordinates=[F(0)]*len(basis)
331 for k in reversed(range(len(basis))):
332 coordinates[k]=remainder[k]
333 for j,value in enumerate(basis[k]):
334 remainder[j]-=coordinates[k]*value
335 assert not any(remainder)
336 return coordinates
337
338
339 def fibonacci_examples(euler_polynomials, signed, check):
340 size=len(euler_polynomials)
341 basis=fibonacci_polynomials(size+1)
342 columns=[monic_coordinates(poly,basis[:size]) for poly in euler_polynomials]
343 connection=[list(row) for row in zip(*columns)]
344 for n in range(size):
345 for k in range(size):
346 closed=sum((euler_polynomials[n][k+2*r]*(-1)**r
347 *(choose(k+2*r,r)-choose(k+2*r,r-1))
348 for r in range(max(0,(n-k)//2+1))),F(0))
349 check('fibonacci_connection_formula',connection[k][n]==closed)
350 recovered=[F(0)]*size
351 for k in range(size):
352 recovered=add(recovered,scale(basis[k],connection[k][n]),size)
353 check('fibonacci_polynomial_reconstruction',trim(recovered)==euler_polynomials[n])
354 recovered_euler=2**n*sum(connection[k][n]*evaluate(basis[k],F(1,2)) for k in range(size))
355 expected=0 if n%2 else signed[n//2]
356 check('fibonacci_centered_euler_recovery',recovered_euler==expected)
357 reverse_columns=[monic_coordinates(poly,euler_polynomials) for poly in basis[:size]]
358 reverse=[list(row) for row in zip(*reverse_columns)]
359 identity=[[F(i==j) for j in range(size)] for i in range(size)]
360 check('fibonacci_connection_inverse',matrix_product(connection,reverse)==identity)
361 check('fibonacci_connection_inverse',matrix_product(reverse,connection)==identity)
362 displayed=[[1,F(-1,2),-1,F(7,4),2],[0,1,-1,-2,5],
363 [0,0,1,F(-3,2),-3],[0,0,0,1,-2],[0,0,0,0,1]]
364 check('fibonacci_display', [row[:5] for row in connection[:5]]==displayed)
365 check('fibonacci_failure_controls',matrix_vector(displayed,[1,0,0,0,0])
366 != matrix_vector(list(map(list,zip(*displayed))),[1,0,0,0,0]))
367 # Direct polynomial matrix iteration versus the Fibonacci recurrence.
368 polynomial_powers=[]
369 state=[[[F(1)],[F(0)]],[[F(0)],[F(1)]]]
370 for n in range(size):
371 polynomial_powers.append(state)
372 if n:
373 expected=[[basis[n],basis[n-1]],
374 [basis[n-1],basis[n-2] if n>=2 else [F(0)]]]
375 check('fibonacci_polynomial_companion',state==expected)
376 state=[[trim(add([F(0)]+state[i][0],state[i][1])),state[i][0]] for i in range(2)]
377
378 count=16
379 a=[F(value) for value in signed[:count]]
380 fib=[0,1]
381 for _ in range(count+1):
382 fib.append(fib[-1]+fib[-2])
383 analysis=[[F(i==j) if i<2 else F(j==i)-F(j==i-1)-F(j==i-2)
384 for j in range(count)] for i in range(count)]
385 synthesis=[]
386 for i in range(count):
387 synthesis.append([F(1 if i==0 else fib[i-1]) if j==0 else
388 F(fib[i]) if j==1 else F(fib[i-j+1]) if j<=i else F(0)
389 for j in range(count)])
390 identity=[[F(i==j) for j in range(count)] for i in range(count)]
391 check('fibonacci_defect_inverse',matrix_product(analysis,synthesis)==identity)
392 check('fibonacci_defect_inverse',matrix_product(synthesis,analysis)==identity)
393 coordinates=matrix_vector(analysis,a)
394 check('fibonacci_defect_reconstruction',matrix_vector(synthesis,coordinates)==a)
395 defects=coordinates[2:]
396 check('fibonacci_display',defects[:4]==[5,-65,1441,-51845])
397 companion=[[F(0),F(1)],[F(1),F(1)]]
398 inverse_companion=[[F(-1),F(1)],[F(1),F(0)]]
399 forward=backward=[[F(1),F(0)],[F(0),F(1)]]
400 moving=[]
401 for n in range(count-1):
402 state=a[n:n+2]
403 z=matrix_vector(backward,state)
404 moving.append(z)
405 check('fibonacci_moving_state',matrix_vector(forward,z)==state)
406 green=matrix_vector(forward,a[:2])
407 for j in range(n):
408 green=add(green,[F(fib[n-1-j]*defects[j]),F(fib[n-j]*defects[j])],2)
409 check('fibonacci_green_reconstruction',green==state)
410 if n<count-2:
411 next_state=add(matrix_vector(companion,state),[0,defects[n]],2)
412 check('fibonacci_forced_companion',next_state==a[n+1:n+3])
413 next_backward=matrix_product(inverse_companion,backward)
414 increment=matrix_vector(next_backward,[0,defects[n]])
415 check('fibonacci_moving_increment',add(z,increment)==matrix_vector(next_backward,next_state))
416 forward=matrix_product(companion,forward)
417 backward=matrix_product(inverse_companion,backward)
418 check('fibonacci_display',moving[:4]==[[1,-1],[6,-1],[71,-66],[2953,-1507]])
419 check('fibonacci_failure_controls',defects[0]!=0)
420
421 residual=mul(a,[F(1),F(-1),F(-1)],count)
422 kernel=inverse(residual,count)
423 transfer=toeplitz_matrix(kernel,count)
424 inverse_transfer=toeplitz_matrix(residual,count)
425 target=[F(fib[n+1]) for n in range(count)]
426 check('fibonacci_sequence_transfer',matrix_vector(transfer,a)==target)
427 check('fibonacci_sequence_transfer',matrix_vector(inverse_transfer,target)==a)
428 check('fibonacci_toeplitz_inverse',matrix_product(transfer,inverse_transfer)==identity)
429 check('fibonacci_toeplitz_inverse',matrix_product(inverse_transfer,transfer)==identity)
430 check('fibonacci_display',residual[:6]==[1,-2,5,-65,1441,-51845])
431 check('fibonacci_display',kernel[:6]==[1,2,-1,53,-1200,46233])
432 # Complete rows of S L(H)-L(H) S: an infinite rank-one boundary term.
433 for n in range(count-1):
434 for j in range(count):
435 first=kernel[n+1-j] if j<=n+1 else 0
436 second=kernel[n-j+1] if 1<=j<=n+1 else 0
437 check('fibonacci_shift_commutator',first-second==(kernel[n+1] if j==0 else 0))
438 return dict(polynomial_connection=[[str(c) for c in row] for row in connection],
439 polynomial_inverse=[[str(c) for c in row] for row in reverse],
440 polynomial_basis=[[str(c) for c in poly] for poly in basis[:size]],
441 polynomial_companion_powers=[[[[str(c) for c in poly] for poly in row] for row in mat] for mat in polynomial_powers],
442 recurrence_analysis=[[str(c) for c in row] for row in analysis],
443 recurrence_synthesis=[[str(c) for c in row] for row in synthesis],
444 defects=[str(c) for c in defects],moving_states=[[str(c) for c in row] for row in moving],
445 ordinary_residual=[str(c) for c in residual],ordinary_kernel=[str(c) for c in kernel],
446 sequence_transfer=[[str(c) for c in row] for row in transfer],
447 sequence_inverse=[[str(c) for c in row] for row in inverse_transfer],
448 target=[int(c) for c in target])
449
450
451 def matrix_power(matrix, exponent):
452 result=[[F(i==j) for j in range(len(matrix))] for i in range(len(matrix))]
453 base=matrix
454 while exponent:
455 if exponent%2:
456 result=matrix_product(result,base)
457 base=matrix_product(base,base)
458 exponent//=2
459 return result
460
461
462 def literature_examples(euler_polynomials, check):
463 fib=[0,1]
464 for _ in range(1400):
465 fib.append(fib[-1]+fib[-2])
466 lucas=[2,1]
467 for _ in range(12):
468 lucas.append(lucas[-1]+lucas[-2])
469 companion=[[F(1),F(1)],[F(1),F(0)]]
470 identity=[[F(1),F(0)],[F(0),F(1)]]
471 jumps=[]
472 for step in range(13):
473 jump=matrix_power(companion,step)
474 jumps.append([[int(c) for c in row] for row in jump])
475 check('literature_lucas_trace',jump[0][0]+jump[1][1]==lucas[step])
476 check('literature_determinant_parity',jump[0][0]*jump[1][1]-jump[0][1]*jump[1][0]==(-1)**step)
477 char=matrix_product(jump,jump)
478 check('literature_jump_characteristic',all(char[i][j]-lucas[step]*jump[i][j]+(-1)**step*identity[i][j]==0
479 for i in range(2) for j in range(2)))
480 for address in range(6):
481 for k in range(9):
482 check('literature_sampled_recurrence',fib[step*(k+2)+address]
483 ==lucas[step]*fib[step*(k+1)+address]-(-1)**step*fib[step*k+address])
484 check('literature_parity_control',fib[5]!=3*fib[3]+fib[1])
485 check('literature_jump_control',matrix_power(companion,2)!=companion)
486
487 precision=[]
488 ladder=[]
489 for refinement in range(5):
490 step=12*2**refinement
491 jump=matrix_power(companion,step)
492 divisor=2**(refinement+3)
493 factor=[[ (jump[i][j]-identity[i][j])/divisor for j in range(2)] for i in range(2)]
494 check('literature_dyadic_refinement',all(c.denominator==1 for row in factor for c in row))
495 check('literature_dyadic_factor_parity',all(int(factor[i][j])%2==int(i==j) for i in range(2) for j in range(2)))
496 for n in range(7):
497 for order in range(6):
498 direct=sum((-1)**(order-j)*comb(order,j)*fib[n+step*j] for j in range(order+1))
499 projected=matrix_vector(matrix_product(matrix_power(companion,n),matrix_power(factor,order)),[1,0])[1]
500 check('literature_matrix_difference',direct==divisor**order*projected)
501 check('literature_dyadic_divisibility',direct%(divisor**order)==0)
502 if n==1:
503 check('literature_dyadic_sharpness',valuation(direct,2)==(refinement+3)*order)
504 for address in range(4):
505 for coarse in range(7):
506 terms=[]
507 for order in range(coarse+1):
508 projected=matrix_vector(matrix_product(matrix_power(companion,address),matrix_power(factor,order)),[1,0])[1]
509 terms.append(comb(coarse,order)*divisor**order*projected)
510 check('literature_packet_binomial',sum(terms)==fib[step*coarse+address])
511 for depth in range(coarse+2):
512 partial=sum(terms[:depth])
513 error=fib[step*coarse+address]-partial
514 check('literature_packet_precision',error%(divisor**depth)==0)
515 if address==1 and coarse==3:
516 precision.append(dict(step=step,refinement=refinement,depth=depth,
517 target=fib[step*coarse+address],partial=int(partial),error=int(error),
518 modulus=divisor**depth,valuation=valuation(error,2)))
519 ladder.append(dict(step=step,power=refinement,divisor=divisor,
520 factor=[[int(c) for c in row] for row in factor]))
521 check('literature_display',ladder[0]['factor']==[[29,18],[18,11]])
522 jump24=matrix_power(companion,24)
523 factor144=[[(jump24[i][j]-identity[i][j])/144 for j in range(2)] for i in range(2)]
524 check('literature_simultaneous_factor',factor144==[[521,322],[322,199]])
525 for n in range(7):
526 for order in range(6):
527 direct=sum((-1)**(order-j)*comb(order,j)*fib[n+24*j] for j in range(order+1))
528 projected=matrix_vector(matrix_product(matrix_power(companion,n),matrix_power(factor144,order)),[1,0])[1]
529 check('literature_simultaneous_differences',direct==144**order*projected)
530 check('literature_simultaneous_moduli',direct%(144**order)==0 and direct%(24**order)==0)
531 check('literature_display',sum((-1)**(2-j)*comb(2,j)*fib[1+24*j] for j in range(3))==7778592000)
532 # Actual composition of the existing ordinary Euler transfer with a selector.
533 a=[(-1)**n*v for n,v in enumerate(reciprocal_recurrence(2,1,16))]
534 kernel=inverse(mul(a,[1,-1,-1],16),16)
535 selector_rows=[toeplitz_matrix(kernel,16)[2*k] for k in range(8)]
536 check('literature_euler_packet_composition',matrix_vector(selector_rows,a)==[fib[2*k+1] for k in range(8)])
537
538 calculus={}
539 for degree in range(9):
540 source=fibonacci_polynomials(degree+1)
541 target=fibonacci_polynomials(degree+2)
542 primitives=[[F(0)]+[value/F(i+1) for i,value in enumerate(poly)] for poly in source]
543 integral=list(map(list,zip(*(monic_coordinates(poly,target) for poly in primitives))))
544 derivative=list(map(list,zip(*(monic_coordinates([F(j)*poly[j] for j in range(1,len(poly))] or [F(0)],source)
545 for poly in target))))
546 c=list(map(list,zip(*(monic_coordinates(poly,source) for poly in euler_polynomials[:degree+1]))))
547 c_inverse=list(map(list,zip(*(monic_coordinates(poly,euler_polynomials[:degree+2]) for poly in target))))
548 d_euler=[[F(j) if i==j-1 else F(0) for j in range(degree+2)] for i in range(degree+1)]
549 conjugate=matrix_product(matrix_product(c,d_euler),c_inverse)
550 check('literature_derivative_conjugacy',conjugate==derivative)
551 combined=matrix_product(conjugate,integral)
552 check('literature_calculus_round_trip',combined==[[F(i==j) for j in range(degree+1)] for i in range(degree+1)])
553 reverse=matrix_product(integral,derivative)
554 evaluation=[poly[0] for poly in target]
555 check('literature_calculus_boundary',reverse==[[F(i==j)-(evaluation[j] if i==0 else 0)
556 for j in range(degree+2)] for i in range(degree+2)])
557 moment=lambda poly:sum(value/F(i+1) for i,value in enumerate(poly))
558 fib_moments=[moment(poly) for poly in source]
559 euler_moments=[moment(poly) for poly in euler_polynomials[:degree+1]]
560 check('literature_integral_covector',matrix_product([fib_moments],c)==[euler_moments])
561 calculus[str(degree)]=dict(integral=integral,derivative=derivative,
562 euler_derivative=d_euler,connection=c,inverse_connection=c_inverse,
563 fib_moments=fib_moments,euler_moments=euler_moments)
564 display=calculus['3']
565 v0=[F(0),F(0),F(0),F(1)]
566 v1=matrix_vector(display['integral'],v0)
567 v2=matrix_vector(display['inverse_connection'],v1)
568 v3=matrix_vector(display['euler_derivative'],v2)
569 v4=matrix_vector(display['connection'],v3)
570 check('literature_display',v1==[F(-1,2),0,F(1,4),0,F(1,4)])
571 check('literature_display',v2==[F(5,8),F(3,2),F(7,4),F(1,2),F(1,4)])
572 check('literature_display',v3==[F(3,2),F(7,2),F(3,2),1] and v4==v0)
573 check('literature_display',sum(x*y for x,y in zip(display['fib_moments'],v4))==F(5,4))
574 check('literature_wrong_operator_control',matrix_vector(display['euler_derivative'],[0,0,0,1,0])
575 !=matrix_vector(display['derivative'],[0,0,0,1,0]))
576 def strings(obj):
577 if isinstance(obj,F):return str(obj)
578 if isinstance(obj,list):return [strings(x) for x in obj]
579 if isinstance(obj,dict):return {k:strings(v) for k,v in obj.items()}
580 return obj
581 return dict(jumps=jumps,lucas=lucas[:13],packet_three=[[fib[3*k+j] for k in range(6)] for j in range(3)],
582 dyadic_factors=ladder,precision=precision,simultaneous_factor=[[int(c) for c in row] for row in factor144],
583 calculus=strings(calculus),calculus_vectors=strings([v0,v1,v2,v3,v4]))
584
585
586 def compose_series(outer, inner, count):
587 result=[F(0)]*count
588 term=[F(1)]+[F(0)]*(count-1)
589 for coefficient in outer[:count]:
590 result=add(result,scale(term,coefficient),count)
591 term=mul(term,inner,count)
592 return result
593
594
595 def stirling_tables(size):
596 first=[[0]*size for _ in range(size)]
597 second=[[0]*size for _ in range(size)]
598 first[0][0]=second[0][0]=1
599 for n in range(1,size):
600 for k in range(1,n+1):
601 first[n][k]=first[n-1][k-1]-(n-1)*first[n-1][k]
602 second[n][k]=second[n-1][k-1]+k*second[n-1][k]
603 return first,second
604
605
606 def cyclic_product(a,b,prime):
607 return [sum(a[j]*b[(i-j)%prime] for j in range(prime))%prime for i in range(prime)]
608
609
610 def atlas_examples(check):
611 gains=[1,1,4,2]
612 original_gains=[1,1,5,3]
613 packets=[]
614 for address in range(4):
615 values=congruential_packet(4,address,48)
616 triangular=[F(1)]
617 for n in range(1,48):
618 triangular.append(-sum(comb(4*n+address,4*k)*triangular[k] for k in range(n))
619 /F(comb(4*n+address,4*n)))
620 check('atlas_quartic_recurrence',triangular==values)
621 check('atlas_quartic_integrality',all(v.denominator%2 for v in values))
622 for n in range(7):
623 for step in range(1,9):
624 for order in range(6):
625 difference=sum((-1)**(order-k)*comb(order,k)*values[n+step*k] for k in range(order+1))
626 weight=valuation(difference,2)
627 check('atlas_all_row_bounds',weight is None or weight>=order*(1+valuation(step,2)))
628 for n in range(9):
629 for step in range(1,17):
630 difference=values[n+step]-values[n]
631 check('atlas_exact_displacement',valuation(difference,2)==valuation(step,2)+gains[address])
632 packets.append(values)
633 newton=[difference_coefficients(values[:13]) for values in packets]
634 for address,orders,expected in [(0,[2],[3]),(1,[2],[3]),(2,[2,3,4,5,6],[6,9,13,16,22]),(3,[2,3,4],[8,8,15])]:
635 check('atlas_newton_prefix',[valuation(newton[address][k],2) for k in orders]==expected)
636 check('atlas_display',[newton[j][1] for j in range(4)]==[F(-2),F(-6,5),F(-16,15),F(-36,35)])
637
638 # All source diagonals are odd, so modular recurrence needs no rational division.
639 modulus=8192
640 binomial=[[comb(4*n,4*k)%modulus for k in range(n+1)] for n in range(256)]
641 modular=[]
642 permutations=[]
643 phi_tables=[]
644 periods=[]
645 for address in range(4):
646 g=[pow(comb(4*k+address,address),-1,modulus) for k in range(256)]
647 values=[1]
648 for n in range(1,256):
649 values.append(-sum(binomial[n][k]*g[k]*values[n-k] for k in range(1,n+1))%modulus)
650 modular.append(values)
651 for n,exact in enumerate(packets[address]):
652 check('atlas_exact_modular_agreement',values[n]==exact.numerator*pow(exact.denominator,-1,modulus)%modulus)
653 factor=factorial(address)
654 d=original_gains[address]
655 original=[factor*v%modulus for v in values]
656 check('atlas_affine_image_residue',all((v-factor)%2**d==0 for v in original))
657 phi=[((v-factor)//2**d)%256 for v in original]
658 phi_tables.append(phi)
659 for bits in range(1,9):
660 precision=2**bits
661 residues=[value%precision for value in phi[:precision]]
662 check('atlas_residue_permutation',sorted(residues)==list(range(precision)))
663 permutations.append(dict(address=address,bits=bits,size=precision,bijection=True))
664 for n in range(32):
665 for m in range(n):
666 check('atlas_finite_isometry',valuation(phi[n]-phi[m],2)==valuation(n-m,2))
667 for bits in range(9):
668 period=2**max(0,bits-d)
669 periods.append(dict(address=address,bits=bits,least_period=period))
670 for n in range(9):
671 check('atlas_least_period',original[n+period]%(2**bits)==original[n]%(2**bits))
672 for step in range(1,17):
673 actual=(original[n+step]-original[n])%(2**bits)==0
674 check('atlas_period_classification',actual==(step%period==0))
675 if period>1:
676 check('atlas_half_period_control',(original[period//2]-original[0])%(2**bits)!=0)
677 trace=[]
678 address,target,index=2,13,0
679 phi=phi_tables[address]
680 for bits in range(6):
681 low=index
682 high=index+2**bits
683 mod=2**(bits+1)
684 candidates=[value for value in (low,high) if phi[value]%mod==target%mod]
685 check('atlas_binary_lifting',len(candidates)==1)
686 index=candidates[0]
687 trace.append(dict(bits=bits+1,modulus=mod,low=low,low_image=phi[low]%mod,
688 high=high,high_image=phi[high]%mod,target=target%mod,chosen=index))
689 check('atlas_inverse_search',[n for n in range(64) if phi[n]%64==13]==[index])
690 row_modulus=2048
691 diagonal=comb(4*index+address,4*index)
692 summands=[comb(4*index+address,4*k)*modular[address][k]%row_modulus for k in range(index)]
693 row_sum=sum(summands)%row_modulus
694 recovered_c=-row_sum*pow(diagonal,-1,row_modulus)%row_modulus
695 recovered_u=factorial(address)*recovered_c%row_modulus
696 check('atlas_inverse_row_witness',recovered_c==modular[address][index]%row_modulus)
697 check('atlas_inverse_row_witness',recovered_u==418)
698 row_witness=dict(index=index,address=address,modulus=row_modulus,diagonal=diagonal,
699 inverse_diagonal=pow(diagonal,-1,row_modulus),summands=summands,
700 sum=row_sum,normalized=recovered_c,original=recovered_u)
701
702 signed_even=[(-1)**n*v for n,v in enumerate(reciprocal_recurrence(2,1,13))]
703 full=[int(signed_even[n//2]) if n%2==0 else 0 for n in range(25)]
704 modes=[]
705 for prime in (3,5,7):
706 unit=[1]+[0]*(prime-1)
707 source=[0]*prime
708 source[1]=source[-1]=pow(2,-1,prime)
709 nilpotent=[(x-y)%prime for x,y in zip(source,unit)]
710 term=unit
711 inverse_weights=[0]*prime
712 for r in range(prime):
713 inverse_weights=[(a+(-1)**r*b)%prime for a,b in zip(inverse_weights,term)]
714 term=cyclic_product(term,nilpotent,prime)
715 check('atlas_mode_nilpotence',term==[0]*prime)
716 check('atlas_mode_inverse',cyclic_product(source,inverse_weights,prime)==unit)
717 check('atlas_mode_inverse',cyclic_product(inverse_weights,source,prime)==unit)
718 values=[sum(weight*pow(u,n,prime) for u,weight in enumerate(inverse_weights))%prime for n in range(25)]
719 check('atlas_modes_euler_values',values==[v%prime for v in full])
720 q=[(1+(-1)**n)*pow(2,-1,prime)%prime for n in range(25)]
721 for n in range(25):
722 product=sum(comb(n,k)*q[k]*values[n-k] for k in range(n+1))%prime
723 check('atlas_modes_hurwitz_product',product==int(n==0))
724 modes.append(dict(prime=prime,source=source,weights=inverse_weights,values=values))
725 check('atlas_display',modes[1]['weights']==[1,1,4,4,1])
726 check('atlas_positive_tail_control',full[0]%5!=full[4]%5)
727 check('atlas_factorial_nonunit',factorial(5)%5==0)
728
729 size=17
730 first,second=stirling_tables(size)
731 identity=[[int(i==j) for j in range(size)] for i in range(size)]
732 check('atlas_stirling_inverse',matrix_product(first,second)==identity)
733 check('atlas_stirling_inverse',matrix_product(second,first)==identity)
734 log=[F(0)]+[F((-1)**(n+1),n) for n in range(1,size)]
735 exponential=[F(0)]+[F(1,factorial(n)) for n in range(1,size)]
736 identity_series=[F(0),F(1)]+[F(0)]*(size-2)
737 check('atlas_log_exp_inverse',compose_series(log,exponential,size)==identity_series)
738 check('atlas_log_exp_inverse',compose_series(exponential,log,size)==identity_series)
739 egf=[F(full[n],factorial(n)) for n in range(size)]
740 logarithmic=compose_series(egf,log,size)
741 transformed=matrix_vector(first,full[:size])
742 check('atlas_log_stirling',transformed==[logarithmic[n]*factorial(n) for n in range(size)])
743 check('atlas_stirling_recovery',matrix_vector(second,transformed)==full[:size])
744 rational=mul([1,1],inverse([1,1,F(1,2)],size),size)
745 check('atlas_log_rational_form',logarithmic==rational)
746 shift=compose_series([F(full[n+1],factorial(n)) for n in range(size)],log,size)
747 extra_log=log+[F((-1)**(size+1),size)]
748 extra_composed=compose_series([F(full[n],factorial(n)) for n in range(size+1)],extra_log,size+1)
749 extra_rational=mul([1,1],inverse([1,1,F(1,2)],size+1),size+1)
750 check('atlas_log_extra_degree',extra_composed[-1]==extra_rational[-1])
751 extended=transformed+[extra_rational[-1]*factorial(size)]
752 theta=[extended[n+1]+n*transformed[n] for n in range(size)]
753 check('atlas_transported_derivative',theta==[shift[n]*factorial(n) for n in range(size)])
754 for n in range(2,size):
755 check('atlas_log_integer_recurrence',transformed[n]==-n*transformed[n-1]-comb(n,2)*transformed[n-2])
756 check('atlas_display',transformed[:9]==[1,0,-1,3,-6,0,90,-630,2520])
757 return dict(quartic=[[str(c) for c in values] for values in packets],newton_prefix=[[str(c) for c in row] for row in newton],
758 gains=gains,original_gains=original_gains,modular_coefficients=modular,modulus=modulus,
759 phi_residues=phi_tables,permutations=permutations,periods=periods,inverse_trace=trace,
760 inverse_address=2,inverse_target=13,inverse_bits=6,inverse_index=index,modes=modes,
761 inverse_row_witness=row_witness,
762 stirling_first=first,stirling_second=second,full_euler=full,logarithmic_euler=[int(v) for v in transformed],
763 logarithmic_ordinary=[str(v) for v in logarithmic],derivative_pullback=[str(v) for v in shift])
764
765
766 def edge_examples(check):
767 kernel=local_kernel(2,70)
768 beta=[kernel[m]/(4**m*factorial(2*m)) for m in range(70)]
769 euler=[(-1)**n*v for n,v in enumerate(reciprocal_recurrence(2,1,70))]
770 carries=[]
771 for r in range(1,256):
772 cutoff=(r.bit_count()-1)//3
773 weights=[3*q+(r+q).bit_count()-r.bit_count() for q in range(cutoff+1)]
774 minimum=min(weights)
775 addresses=[q for q,w in enumerate(weights) if w==minimum]
776 direct=[3*q+(r+q).bit_count()-r.bit_count() for q in range(257)]
777 check('edge_carry_cutoff',minimum==min(direct) and addresses==[q for q,w in enumerate(direct) if w==minimum])
778 if r in (1,3,7,15,31,63,127,254,255):
779 reduction=[0]
780 for q in addresses:
781 polynomial=[1]
782 for i in range(q):polynomial=mul(polynomial,[-i,1])
783 reduction=add(reduction,polynomial)
784 carries.append(dict(r=r,binary_ones=r.bit_count(),cutoff=cutoff,mu=minimum,
785 minimizing_addresses=addresses,reduction_mod_2=[int(c)%2 for c in trim(reduction)]))
786 rows=[]
787 for r,offsets in [(r,range(9)) for r in range(1,9)]+[(r,range(7)) for r in (15,31,63)]:
788 for x in offsets:
789 n=r+x
790 terms=[kernel[m]*comb(2*n,2*m)*local_source(2,2*(n-m)) for m in range(n+1)]
791 check('edge_direct_euler_reconstruction',sum(terms)==euler[n])
792 partial=sum(terms[:r])
793 tail=euler[n]-partial
794 h=F(0)
795 polynomial=1
796 for q in range(x+1):
797 if q:polynomial*=(x-q+1)*(2*x-(2*q-1))
798 h+=8**q*beta[r+q]/beta[r]*polynomial*local_source(2,2*(x-q))
799 factor=kernel[r]*comb(2*n,2*r)
800 check('edge_finite_factorization',tail==factor*h)
801 hv=valuation(h,2)
802 tv=valuation(tail,2)
803 binomial_weight=valuation(comb(n,r),2)
804 check('edge_binomial_weight',binomial_weight==valuation(comb(2*n,2*r),2))
805 check('edge_valuation_identity',tv is None if hv is None else tv==4*r-1+binomial_weight+hv)
806 if r.bit_count()<=3:
807 check('edge_sparse_binary_family',hv==0 and tv==4*r-1+binomial_weight)
808 if x==0:check('edge_uniform_sharpness',h==1 and tail==kernel[r] and tv==4*r-1)
809 rows.append(dict(r=r,N=n,x=x,retained_depth=r-1,partial=str(partial),tail=str(tail),
810 H=str(h),H_valuation=hv,tail_valuation=tv,binomial_valuation=binomial_weight,
811 minimum_summand_valuation=min(valuation(t,2) for t in terms[r:]),
812 omitted_weights=[valuation(t,2) for t in terms[r:]]))
813 witness=next(row for row in rows if row['r']==15 and row['x']==1)
814 check('edge_cancellation_witness',witness['tail_valuation']>witness['minimum_summand_valuation'])
815 return dict(carries=carries,rows=rows,cancellation_witness=witness,
816 kernel={str(r):str(kernel[r]) for r in list(range(1,9))+[15,31]})
817
818
819 def compute():
820 counts = {}
821
822 def check(group, condition):
823 if not condition:
824 raise AssertionError('Failed: '+group)
825 counts[group] = counts.get(group,0)+1
826
827 # Route A: compressed EGF recurrence versus ordinary rational products.
828 parameters = [F(-2), F(-1), F(0), F(1,2), F(1), F(2), F(3)]
829 for r in parameters:
830 for n,value in enumerate(reciprocal_recurrence(1,r,11)):
831 check('spacing_one_boundary',value == r**n)
832 for spacing in range(2,7):
833 polynomials = reciprocal_polynomials(spacing,11)
834 for r in parameters:
835 a = reciprocal_recurrence(spacing,r,11)
836 b = reciprocal_binomial(spacing,r,11)
837 for n in range(11):
838 check('reciprocal_routes',a[n] == b[n])
839 check('parameter_polynomials',evaluate(polynomials[n],r) == a[n])
840 # Compare every coefficient of the two-variable addition identity.
841 for n in range(9):
842 left, right = {}, {}
843 for degree,c in enumerate(polynomials[n]):
844 for j in range(degree+1):
845 left[j,degree-j] = left.get((j,degree-j),0)+c*comb(degree,j)
846 for k in range(n+1):
847 for i,a in enumerate(polynomials[k]):
848 for j,b in enumerate(polynomials[n-k]):
849 right[i,j] = right.get((i,j),0)+comb(spacing*n,spacing*k)*a*b
850 clean = lambda x:{k:v for k,v in x.items() if v}
851 check('addition_polynomial_identities',clean(left) == clean(right))
852
853 s = reciprocal_recurrence(2,1,41)
854 b = reciprocal_recurrence(4,1,21)
855 for n,value in enumerate(SECANT):
856 check('oeis_secant_prefix',s[n] == value)
857 for n,value in enumerate(QUARTIC):
858 check('oeis_quartic_prefix',b[n] == value)
859 for spacing in range(2,7):
860 a = reciprocal_recurrence(spacing,1,9)
861 for n in range(9):
862 check('rank_dynamic_program',a[n] == alternating_count(spacing,spacing*n))
863 check('permutation_composition_formula',a[n] == composition_count(spacing,n))
864 # An explicit enumeration checks the rank-DP implementation on small inputs.
865 for spacing,size in [(2,0),(2,4),(2,6),(2,8),(4,4),(4,8)]:
866 actual = sum(all((p[j-1] < p[j]) == (j % spacing == 0)
867 for j in range(1,size))
868 for p in permutations(range(size)))
869 check('permutation_enumeration',actual == alternating_count(spacing,size))
870
871 # Route B: the 0/2 mod 4 transfer matrix and its determinant.
872 length = 12
873 c = [F(1,factorial(4*n)) for n in range(length)]
874 d = [-F(1,factorial(4*n+2)) for n in range(length)]
875 g = add(mul(c,c,length),scale(shift(mul(d,d,length),1,length),-1),length)
876 expected = [F((-4)**n,factorial(4*n)) for n in range(length)]
877 for x,y in zip(g,expected):
878 check('quartic_determinant',x == y)
879 h = inverse(g,length)
880 packet_a = mul(c,h,length)
881 packet_b = scale(mul(d,h,length),-1)
882 first = add(mul(c,packet_a,length),shift(mul(d,packet_b,length),1,length),length)
883 second = add(mul(d,packet_a,length),mul(c,packet_b,length),length)
884 for n in range(length):
885 check('packet_reconstruction',packet_a[n]*factorial(4*n) == s[2*n])
886 check('packet_reconstruction',packet_b[n]*factorial(4*n+2) == s[2*n+1])
887 check('packet_product_certificate',first[n] == (n == 0) and second[n] == 0)
888 for n in range(41):
889 rebuilt = sum(comb(2*n,4*k)*4**k*b[k] for k in range(n//2+1))
890 check('quartic_to_secant',rebuilt == s[n])
891 recovered=[]
892 for n in range(21):
893 remainder=s[2*n]-sum(comb(4*n,4*k)*4**k*recovered[k] for k in range(n))
894 recovered.append(remainder/4**n)
895 check('secant_to_quartic',recovered[n] == b[n])
896 naive_a = inverse(c,length)
897 check('failure_controls',naive_a[1] == -F(1,24) and packet_a[1] == F(5,24))
898 check('failure_controls',mul(d,naive_a,length)[0] == -F(1,2))
899 check('failure_controls',packet_a[1]-naive_a[1] == F(1,4))
900 check('failure_controls',F(s[2],factorial(2)) != F(s[2],factorial(4)))
901 check('failure_controls',add(mul(c,c,2),scale(mul(d,d,2),-1),2)[0] != g[0])
902
903 # Route C: Newton residual squaring in Q[t]/(t^32).
904 precision = 32
905 cosine = [F(0)]*precision
906 for n in range(precision//2):
907 cosine[2*n] = F((-1)**n,factorial(2*n))
908 exact_inverse = inverse(cosine,precision)
909 for n in range(precision//2):
910 check('dense_scalar_inverse',exact_inverse[2*n]*factorial(2*n) == s[n])
911 current = [F(1)]+[F(0)]*(precision-1)
912 residual_orders = []
913 for iteration in range(5):
914 residual = add([1],scale(mul(cosine,current,precision),-1),precision)
915 first_nonzero = next((i for i,value in enumerate(residual) if value),None)
916 expected_order = min(2**(iteration+1),precision)
917 check('newton_precision',first_nonzero == (expected_order if expected_order < precision else None))
918 residual_orders.append(first_nonzero if first_nonzero is not None else 'at least 32')
919 check('newton_precision',current[:expected_order] == exact_inverse[:expected_order])
920 next_value = mul(current,add([2],scale(mul(cosine,current,precision),-1),precision),precision)
921 next_residual = add([1],scale(mul(cosine,next_value,precision),-1),precision)
922 check('newton_residual_identity',next_residual == mul(residual,residual,precision))
923 current = next_value
924 finite_s3 = [F(0)]*10
925 for n in range(4):
926 finite_s3[2*n] = s[n]/factorial(2*n)
927 actual_defect = add([1],scale(mul(cosine,finite_s3,10),-1),10)
928 check('truncation_defect',actual_defect[:8] == [0]*8 and actual_defect[8] == F(277,8064))
929
930 # Euler polynomials use a different variable and the signed sech convention.
931 size = 13
932 denominator = [F(2)]+[F(1,factorial(n)) for n in range(1,size)]
933 base = scale(inverse(denominator,size),2)
934 euler_polynomials = []
935 signed_full = [F(0)]*size
936 for n in range(0,size,2):
937 signed_full[n] = (-1)**(n//2)*s[n//2]
938 for n in range(size):
939 poly = [comb(n,j)*base[n-j]*factorial(n-j) for j in range(n+1)]
940 euler_polynomials.append(trim(poly))
941 centered = [F(0)]*(n+1)
942 for k in range(n//2+1):
943 piece = translate([F(0)]*(n-2*k)+[F(1)],-F(1,2))
944 centered = add(centered,scale(piece,comb(n,2*k)*signed_full[2*k]/2**(2*k)),n+1)
945 check('euler_polynomial_centering',trim(poly) == trim(centered))
946 target = [F(0)]*n+[F(2)]
947 check('euler_polynomial_shift',trim(add(translate(poly,1),poly)) == target)
948 check('euler_polynomial_reflection',translate(poly,1,-1) == trim(scale(poly,(-1)**n)))
949 if n:
950 derivative = [j*poly[j] for j in range(1,len(poly))]
951 check('euler_polynomial_derivative',trim(derivative) == trim(scale(euler_polynomials[n-1],n)))
952 for length in (0,1,2,3,5,6):
953 for start in (-F(1,2),F(0),F(1),F(3,2)):
954 direct=sum((-1)**j*(start+j)**n for j in range(length))
955 boundary=(evaluate(poly,start)-(-1)**length*evaluate(poly,start+length))/2
956 check('alternating_power_telescope',direct == boundary)
957
958 # Paper IX charts: direct windows and rational sources versus formal units.
959 local_count=11
960 bernoulli=bernoulli_numbers(20)
961 signed_egf=[(-1)**n*s[n]/factorial(2*n) for n in range(local_count)]
962 for chart in (1,2,3,5,7,9):
963 kernel=local_kernel(chart,local_count)
964 inverse_coefficients=local_inverse_kernel(chart,local_count)
965 ordinary_inverse=[v/factorial(2*r) for r,v in enumerate(inverse_coefficients)]
966 source_by_series=mul(ordinary_inverse,signed_egf,local_count)
967 for n in range(local_count):
968 source=local_source(chart,2*n)
969 check('local_source_series',source_by_series[n]*factorial(2*n) == source)
970 rebuilt=sum(local_terms(chart,n))
971 check('local_euler_reconstruction',rebuilt == (-1)**n*s[n])
972 inverse_source=sum(comb(2*n,2*r)*inverse_coefficients[r]*(-1)**(n-r)*s[n-r]
973 for r in range(n+1))
974 check('local_inverse_reconstruction',inverse_source == source)
975 if chart == 2:
976 check('dyadic_source_units',valuation(source,2) == 0)
977 if n:
978 check('dyadic_bernoulli_kernel',kernel[n] == -2**(4*n)*(2*n-1)*bernoulli[2*n])
979 check('dyadic_kernel_valuation',valuation(kernel[n],2) == 4*n-1)
980 else:
981 check('odd_reciprocal_power_kernel',kernel[n] == (-1)**n*chart**(2*n)*reciprocal_recurrence(2,2,n+1)[n])
982 terms=local_terms(chart,n)
983 for depth in range(n+1):
984 error=rebuilt-sum(terms[:depth+1])
985 if chart in (2,3,5,7):
986 bound=4*depth+3 if chart == 2 else 2*depth+2
987 actual=valuation(error,chart)
988 check('local_truncation_bound',actual is None or actual >= bound)
989 if depth == n:
990 check('full_depth_exact',error == 0)
991 if chart == 2:
992 for depth in range(10):
993 error=sum(local_terms(2,depth+1))-sum(local_terms(2,depth+1)[:depth+1])
994 check('dyadic_first_omitted_sharpness',error == kernel[depth+1] and valuation(error,2) == 4*depth+3)
995 odd_example=local_row_record(3,4)
996 dyadic_example=local_row_record(2,4)
997 check('displayed_local_reconstruction',odd_example['exact'] == dyadic_example['exact'] == '1385')
998 check('displayed_local_reconstruction',[r['valuation'] for r in dyadic_example['rows']] == [5,8,13,15,None])
999 check('displayed_local_reconstruction',[r['term'] for r in dyadic_example['rows']]
1000 == ['7381/45','-6560/3','163072/15','-204800/9','229376/15'])
1001 check('displayed_local_reconstruction',rational_residue(F(7381,45),32) == 1385 % 32 == 9)
1002 try:
1003 rational_residue(F(5,24),8)
1004 except ValueError:
1005 check('failure_controls',True)
1006 else:
1007 check('failure_controls',False)
1008 check('failure_controls',sum((-1)**a for a in range(2)) == 0)
1009
1010 # Infinite transforms tested through their declared, exact finite blocks.
1011 matrix_size=12
1012 identity=[[F(n == j) for j in range(matrix_size)] for n in range(matrix_size)]
1013 matrices={}
1014 for chart in (2,3,5,7):
1015 kernel=local_kernel(chart,matrix_size)
1016 inverse_coefficients=local_inverse_kernel(chart,matrix_size)
1017 forward=coefficient_matrix(kernel,matrix_size)
1018 backward=coefficient_matrix(inverse_coefficients,matrix_size)
1019 matrices[chart]=forward
1020 for product in (matrix_product(forward,backward),matrix_product(backward,forward)):
1021 for n in range(matrix_size):
1022 for j in range(matrix_size):
1023 check('local_matrix_inverse_entries',product[n][j] == identity[n][j])
1024 source=[local_source(chart,2*n) for n in range(matrix_size)]
1025 target=[(-1)**n*s[n] for n in range(matrix_size)]
1026 for actual,expected_value in zip(matrix_vector(forward,source),target):
1027 check('local_matrix_reconstruction',actual == expected_value)
1028 for actual,expected_value in zip(matrix_vector(backward,target),source):
1029 check('local_matrix_inverse_reconstruction',actual == expected_value)
1030 ordinary=toeplitz_matrix([v/factorial(2*n) for n,v in enumerate(kernel)],matrix_size)
1031 for n in range(matrix_size):
1032 for j in range(matrix_size):
1033 check('factorial_matrix_conjugation',forward[n][j]
1034 == F(factorial(2*n),factorial(2*j))*ordinary[n][j])
1035 for length in range(1,matrix_size+1):
1036 block=lambda a:[row[:length] for row in a[:length]]
1037 check('stable_triangular_prefixes',matrix_product(block(forward),block(backward)) == block(identity))
1038 for left,right in ((2,3),(3,5),(5,2)):
1039 a,b_kernel=local_kernel(left,matrix_size),local_kernel(right,matrix_size)
1040 convolution=[sum(comb(2*n,2*r)*a[r]*b_kernel[n-r] for r in range(n+1))
1041 for n in range(matrix_size)]
1042 product=matrix_product(matrices[left],matrices[right])
1043 combined=coefficient_matrix(convolution,matrix_size)
1044 for n in range(matrix_size):
1045 for j in range(matrix_size):
1046 check('coefficient_matrix_composition',product[n][j] == combined[n][j])
1047 p,q=pascal_matrix(matrix_size),pascal_matrix(matrix_size,inverse=True)
1048 product=matrix_product(p,q)
1049 for n in range(matrix_size):
1050 for j in range(matrix_size):
1051 check('pascal_inverse_entries',product[n][j] == identity[n][j])
1052 euler_vector=[as_integer((-1)**n*s[n]) for n in range(matrix_size)]
1053 differences=difference_coefficients(euler_vector)
1054 for actual,expected_value in zip(matrix_vector(q,euler_vector),differences):
1055 check('pascal_coefficient_coordinates',actual == expected_value)
1056 for actual,expected_value in zip(matrix_vector(p,differences),euler_vector):
1057 check('pascal_coefficient_coordinates',actual == expected_value)
1058 for spacing in range(1,4):
1059 for offset in range(spacing):
1060 count=3*spacing+offset+1
1061 selector=[[F(j == spacing*k+offset) for j in range(count)] for k in range(4)]
1062 conjugated=matrix_product(matrix_product(pascal_matrix(4,inverse=True),selector),pascal_matrix(count))
1063 for k in range(4):
1064 weights=transfer_row(spacing,offset,k)
1065 expected_row=list(weights)+[F(0)]*(count-len(weights))
1066 for actual,expected_value in zip(conjugated[k],expected_row):
1067 check('sampling_pascal_conjugation',actual == expected_value)
1068 lane=[euler_vector[spacing*k+offset] for k in range(4)]
1069 direct=difference_coefficients(lane)
1070 for actual,expected_value in zip(matrix_vector(conjugated,differences[:count]),direct):
1071 check('sampling_matrix_application',actual == expected_value)
1072 sampling=[list(transfer_row(2,0,k))+[F(0)]*(5-len(transfer_row(2,0,k))) for k in range(3)]
1073 correct_cutoff=matrix_vector(sampling,differences[:5])
1074 wrong_cutoff=matrix_vector([row[:3] for row in sampling],differences[:3])
1075 check('sampling_cutoff_counterexample',correct_cutoff == [1,4,1376] and wrong_cutoff == [1,4,32])
1076 wrong_ordinary=toeplitz_matrix([v/factorial(2*r) for r,v in enumerate(local_kernel(3,2))],2)
1077 check('matrix_normalization_counterexample',matrix_vector(wrong_ordinary,[1,17]) == [1,8])
1078
1079 # Regroup a scalar infinite Toeplitz operator into two packet blocks.
1080 packet_size=6
1081 cosine_coefficients=[F((-1)**n,factorial(2*n)) for n in range(2*packet_size)]
1082 scalar_operator=toeplitz_matrix(cosine_coefficients,2*packet_size)
1083 permutation=list(range(0,2*packet_size,2))+list(range(1,2*packet_size,2))
1084 regrouped=[[scalar_operator[n][j] for j in permutation] for n in permutation]
1085 packet_c=[F(1,factorial(4*n)) for n in range(packet_size)]
1086 packet_d=[-F(1,factorial(4*n+2)) for n in range(packet_size)]
1087 cc=toeplitz_matrix(packet_c,packet_size)
1088 dd=toeplitz_matrix(packet_d,packet_size)
1089 ud=toeplitz_matrix([F(0)]+packet_d[:-1],packet_size)
1090 block_operator=[cc[n]+ud[n] for n in range(packet_size)]+[dd[n]+cc[n] for n in range(packet_size)]
1091 for n in range(2*packet_size):
1092 for j in range(2*packet_size):
1093 check('infinite_packet_block_entries',regrouped[n][j] == block_operator[n][j])
1094 original=[s[n]/factorial(2*n) for n in range(2*packet_size)]
1095 packet_vector=[original[n] for n in permutation]
1096 for n,value in enumerate(matrix_vector(block_operator,packet_vector)):
1097 check('packet_block_inverse_application',value == (n == 0))
1098 matrix_examples=dict(local_size=matrix_size,
1099 odd_three=[[str(v) for v in row[:5]] for row in matrices[3][:5]],
1100 dyadic=[[str(v) for v in row[:5]] for row in matrices[2][:5]],
1101 odd_three_inverse=[[str(v) for v in row[:5]] for row in coefficient_matrix(local_inverse_kernel(3,matrix_size),matrix_size)[:5]],
1102 sampling_two=[[str(v) for v in row] for row in sampling],
1103 correct_sampling_cutoff=[str(v) for v in correct_cutoff],wrong_sampling_cutoff=[str(v) for v in wrong_cutoff])
1104
1105 # Exact base changes between local source coordinates, before modular reduction.
1106 base_size=8
1107 charts=(2,3,5,7)
1108 base_changes={}
1109 for source_chart in charts:
1110 for target_chart in charts:
1111 inverse_target=coefficient_matrix(local_inverse_kernel(target_chart,base_size),base_size)
1112 forward_source=[row[:base_size] for row in matrices[source_chart][:base_size]]
1113 change=matrix_product(inverse_target,forward_source)
1114 base_changes[source_chart,target_chart]=change
1115 source_vector=[local_source(source_chart,2*n) for n in range(base_size)]
1116 target_vector=[local_source(target_chart,2*n) for n in range(base_size)]
1117 for actual,expected_value in zip(matrix_vector(change,source_vector),target_vector):
1118 check('local_base_change_values',actual == expected_value)
1119 for row in change:
1120 for value in row:
1121 check('dyadic_integral_base_change',value.denominator % 2 == 1)
1122 if source_chart != 2 and target_chart != 2:
1123 check('odd_integral_base_change',value.denominator == 1)
1124 for source_chart in charts:
1125 for intermediate in charts:
1126 for target_chart in charts:
1127 composed=matrix_product(base_changes[intermediate,target_chart],base_changes[source_chart,intermediate])
1128 direct=base_changes[source_chart,target_chart]
1129 for n in range(base_size):
1130 for j in range(base_size):
1131 check('base_change_cocycle_entries',composed[n][j] == direct[n][j])
1132 loss_entry=base_changes[3,2][4][2]
1133 perturbation=[F(0)]*base_size
1134 perturbation[2]=9
1135 output_error=matrix_vector(base_changes[3,2],perturbation)[4]
1136 check('base_change_precision_loss',loss_entry == F(215264,3))
1137 check('base_change_precision_loss',output_error == 645792 and valuation(output_error,3) == 1)
1138 check('base_change_precision_loss',valuation(loss_entry,3) == -1 and valuation(loss_entry,2) == 5)
1139 ladders=[]
1140 for chart in (2,3):
1141 rows=local_row_record(chart,4)['rows']
1142 for depth in range(4):
1143 modulus=chart**(4*depth+3 if chart==2 else 2*depth+2)
1144 partial=F(rows[depth]['partial'])
1145 residue=rational_residue(partial,modulus)
1146 check('displayed_supercongruence_ladders',residue == 1385 % modulus)
1147 ladders.append(dict(chart=chart,depth=depth,modulus=modulus,partial=str(partial),
1148 residue=residue,error=rows[depth]['error'],valuation=rows[depth]['valuation']))
1149 odd_source=[local_source(3,2*n) for n in range(5)]
1150 changed=matrix_vector([row[:5] for row in base_changes[3,5][:5]],odd_source)
1151 check('odd_chart_modulus_transport',rational_residue(changed[2],81) == 7)
1152 check('odd_chart_modulus_transport',F(896)+192*17+545 == changed[2] == 4705)
1153 check('modulus_change_counterexample',1 % 81 == 82 % 81 and 1 % 25 != 82 % 25)
1154 crt_residue=9+32*((8-9)*pow(32,-1,81) % 81)
1155 check('coprime_supercongruence_combination',32*38-81*15 == 1)
1156 check('coprime_supercongruence_combination',crt_residue == 1385 and crt_residue % 81 == 8 and crt_residue % 32 == 9)
1157 for spacing in range(1,7):
1158 for address in range(spacing):
1159 packet=congruential_packet(spacing,address,9)
1160 triangular=[F(1)]
1161 clearing=1
1162 for n in range(9):
1163 if n:
1164 diagonal=comb(spacing*n+address,spacing*n)
1165 triangular.append(-sum(comb(spacing*n+address,spacing*j)*triangular[j] for j in range(n))/F(diagonal))
1166 clearing*=diagonal
1167 check('congruential_family_routes',packet[n] == triangular[n])
1168 check('congruential_denominator_clearing',(clearing*packet[n]).denominator == 1)
1169 family_four=congruential_packet(4,0,9)
1170 family_two=congruential_packet(2,0,9)
1171 family_boundary=congruential_packet(2,1,9)
1172 for n in range(9):
1173 check('congruential_family_specializations',family_four[n] == (-1)**n*b[n])
1174 check('congruential_family_specializations',family_two[n] == (-1)**n*s[n])
1175 check('congruential_family_specializations',family_boundary[n] == (2-2**(2*n))*bernoulli[2*n])
1176 base_change_examples=dict(size=base_size,
1177 three_to_five=[[str(v) for v in row[:5]] for row in base_changes[3,5][:5]],
1178 three_to_two=[[str(v) for v in row[:5]] for row in base_changes[3,2][:5]],
1179 three_adic_loss_entry=str(loss_entry),input_error=9,output_error=str(output_error),
1180 congruence_ladders=ladders,crt=dict(moduli=[81,32],residues=[8,9],residue=crt_residue,modulus=81*32),
1181 family_four=[str(v) for v in family_four[:5]])
1182
1183 # The same signed Euler data, now in integer finite-difference coordinates.
1184 signed = [as_integer((-1)**n*s[n]) for n in range(41)]
1185 signals = [signed, [n**3-2*n*n+7 for n in range(41)],
1186 [(-1)**n*(n+1)+n*n for n in range(41)]]
1187 for signal in signals:
1188 differences = difference_coefficients(signal)
1189 for spacing in range(1,6):
1190 for offset in range(spacing):
1191 lane = [signal[spacing*q+offset] for q in range(6)]
1192 direct = difference_coefficients(lane)
1193 for order in range(6):
1194 weights = transfer_row(spacing,offset,order)
1195 via_kernel = sum(weight*differences[j] for j,weight in enumerate(weights))
1196 check('finite_difference_transfer',direct[order] == via_kernel)
1197 for m in range(1,5):
1198 for ell in range(1,5):
1199 for offset in range(m):
1200 for branch in range(ell):
1201 for order in range(4):
1202 outer = transfer_row(ell,branch,order)
1203 composed = [F(0)]
1204 for j,weight in enumerate(outer):
1205 composed = add(composed,scale(transfer_row(m,offset,j),weight))
1206 direct = transfer_row(m*ell,offset+m*branch,order)
1207 check('refinement_composition',trim(composed) == direct)
1208 for prime in (2,3,5,7,11):
1209 residual = [choose(prime,j)-(1 if j in (0,prime) else 0) for j in range(prime+1)]
1210 check('prime_residual_polynomials',all(value % prime == 0 for value in residual))
1211 correction = [value//prime for value in residual]
1212 for signal in signals:
1213 differences = difference_coefficients(signal[:prime+1])
1214 check('prime_residual_application',signal[prime]-signal[0]
1215 == differences[prime]+prime*sum(correction[j]*differences[j] for j in range(prime+1)))
1216 c0 = difference_coefficients(signed[:5])
1217 check('displayed_integer_calculations',c0 == [1,-2,8,-80,1664])
1218 check('displayed_integer_calculations',2*c0[1]+c0[2] == 4)
1219 check('displayed_integer_calculations',4*c0[2]+4*c0[3]+c0[4] == 1376)
1220 check('displayed_integer_calculations',4*c0[1]+6*c0[2]+4*c0[3]+c0[4] == 1384)
1221 check('displayed_integer_calculations',2*4+1376 == 1384)
1222 check('failure_controls',c0[2] != signed[2]-signed[0])
1223 check('failure_controls',(signed[2]-signed[0]) % 8 != c0[2] % 8)
1224 check('failure_controls',comb(4,2) % 4 != 0)
1225 first_six=difference_coefficients(signed[:6])
1226 check('higher_prime_examples',first_six == [1,-2,8,-80,1664,-58112])
1227 for prime,expected_correction in ((3,18),(5,7590)):
1228 correction=sum(comb(prime,j)*first_six[j] for j in range(1,prime))
1229 check('higher_prime_examples',correction == expected_correction)
1230 check('higher_prime_examples',first_six[prime]+correction == signed[prime]-1)
1231 alternating=[(-1)**n for n in range(12)]
1232 constant=[1]*12
1233 check('sampled_signal_control',alternating[::2] == constant[::2])
1234 check('sampled_signal_control',alternating[1::2] != constant[1::2])
1235 check('sampled_signal_control',all(alternating[n+2]-alternating[n] == 0 for n in range(10)))
1236 # No field-valued quotient or division by a non-unit is used in the modular check.
1237
1238 fibonacci=fibonacci_examples(euler_polynomials,signed,check)
1239 literature=literature_examples(euler_polynomials,check)
1240 atlas=atlas_examples(check)
1241 edge=edge_examples(check)
1242 polynomials = reciprocal_polynomials(2,6)
1243 report = dict(status='PASS',arithmetic='exact integers and fractions',
1244 python=platform.python_version(),randomness=False,checks=counts,
1245 source_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
1246 secant=[as_integer(v) for v in s[:16]],quartic=[as_integer(v) for v in b[:12]],
1247 signed_euler=signed[:9],difference_coordinates=c0,
1248 lane_spacing_two=[signed[0],signed[2],signed[4]],
1249 lane_difference_coordinates=[1,4,1376],refined_first_difference=1384,
1250 newton_residual_orders=residual_orders,
1251 reciprocal_power_polynomials=[[str(c) for c in p] for p in polynomials],
1252 euler_polynomials=[[str(c) for c in p] for p in euler_polynomials[:7]],
1253 local_reconstructions=[odd_example,dyadic_example],
1254 matrix_examples=matrix_examples,
1255 base_change_examples=base_change_examples,
1256 fibonacci=fibonacci,
1257 literature=literature,
1258 atlas=atlas,
1259 edge=edge,
1260 references=REFERENCES,
1261 bounds=dict(reciprocal_spacing=[2,6],reciprocal_index=[0,10],
1262 reciprocal_parameters=[str(r) for r in parameters],
1263 addition_spacing=[2,6],addition_index=[0,8],
1264 quartic_bridge_secant_index=[0,40],quartic_packet_index=[0,11],
1265 permutation_dp_spacing=[2,6],permutation_dp_block_count=[0,8],
1266 enumeration_cases=[[2,0],[2,4],[2,6],[2,8],[4,4],[4,8]],
1267 euler_polynomial_degree=[0,12],newton_precision=32,
1268 alternating_power_lengths=[0,1,2,3,5,6],alternating_power_starts=['-1/2','0','1','3/2'],
1269 local_charts=[1,2,3,5,7,9],local_even_degree=[0,20],
1270 dyadic_sharp_depth=[0,9],bernoulli_degree=[0,20],
1271 local_matrix_size=12,packet_block_size=6,pascal_matrix_size=12,
1272 sampling_matrix_spacings=[1,3],sampling_matrix_orders=[0,3],
1273 base_change_charts=[2,3,5,7],base_change_matrix_size=8,base_change_all_triples=True,
1274 congruential_family_spacing=[1,6],congruential_addresses='every 0<=j<N',congruential_index=[0,8],
1275 fibonacci_polynomial_degree=[0,12],fibonacci_sequence_length=16,
1276 fibonacci_companion_power=[0,12],fibonacci_state_index=[0,14],
1277 literature_packet_step=[0,12],literature_packet_address=[0,5],literature_packet_index=[0,8],
1278 literature_dyadic_refinement=[0,4],literature_difference_base=[0,6],literature_difference_order=[0,5],
1279 literature_binomial_address=[0,3],literature_binomial_index=[0,6],literature_calculus_degree=[0,8],
1280 quartic_exact_index=[0,47],quartic_difference_base=[0,6],quartic_difference_step=[1,8],quartic_difference_order=[0,5],
1281 quartic_displacement_base=[0,8],quartic_displacement_step=[1,16],quartic_modular_index=[0,255],quartic_modulus=8192,
1282 residue_permutation_bits=[1,8],isometry_pair_indices=[0,31],period_precision=[0,8],inverse_bits=6,
1283 mode_primes=[3,5,7],mode_full_degree=[0,24],stirling_degree=16,log_derivative_input_degree=17,
1284 edge_small_r=[1,8],edge_small_x=[0,8],edge_large_r=[15,31,63],edge_large_x=[0,6],edge_max_full_degree=138,
1285 edge_carry_r=[1,255],edge_carry_q=[0,256],
1286 difference_spacing=[1,5],difference_offsets='all 0 <= s < M',difference_order=[0,5],
1287 composition_spacings=[1,4],composition_offsets='all admitted offsets',composition_order=[0,3],
1288 primes=[2,3,5,7,11]),
1289 non_claims=['No novelty or exclusive-method claim.',
1290 'Finite computations are separate from the written universal algebraic proofs.',
1291 'No numerical error bound, asymptotic law or complexity advantage is established.',
1292 'No Lean execution. The local truncation theorems are proved separately; finite checks do not certify broader Euler congruences.'])
1293 return report
1294
1295
1296 def save_results(directory, report):
1297 directory.mkdir(parents=True,exist_ok=True)
1298 paths = [directory/'results.json',directory/'euler-coefficients.csv',directory/'power-polynomials.csv',
1299 directory/'local-reconstruction.csv',directory/'coefficient-matrices.csv',directory/'congruence-ladders.csv',
1300 directory/'fibonacci-transfer.csv',directory/'fibonacci-precision.csv',directory/'polynomial-calculus.csv',
1301 directory/'quartic-addresses.csv',directory/'dyadic-inverse.csv',directory/'euler-coordinate-atlas.csv',
1302 directory/'dyadic-edge-tails.csv']
1303 if any(p.exists() for p in paths):
1304 raise FileExistsError('Choose a directory with no existing result files')
1305 paths[0].write_text(json.dumps(report,indent=2)+'\n',encoding='utf-8',newline='\n')
1306 with paths[1].open('w',encoding='utf-8',newline='') as stream:
1307 writer = csv.writer(stream,lineterminator='\n')
1308 writer.writerow(['compressed_index','ordinary_degree','secant','signed_euler','ordinary_secant_coefficient'])
1309 for n,s in enumerate(report['secant']):
1310 writer.writerow([n,2*n,s,(-1)**n*s,str(F(s,factorial(2*n)))])
1311 with paths[2].open('w',encoding='utf-8',newline='') as stream:
1312 writer = csv.writer(stream,lineterminator='\n')
1313 writer.writerow(['spacing','compressed_index','parameter_power','coefficient'])
1314 for spacing in (2,4):
1315 for n,p in enumerate(reciprocal_polynomials(spacing,9)):
1316 for degree,value in enumerate(p):
1317 writer.writerow([spacing,n,degree,str(value)])
1318 with paths[3].open('w',encoding='utf-8',newline='') as stream:
1319 writer=csv.writer(stream,lineterminator='\n')
1320 writer.writerow(['chart','ordinary_degree','depth','term','partial','error','error_valuation'])
1321 for example in report['local_reconstructions']:
1322 for row in example['rows']:
1323 writer.writerow([example['chart'],example['ordinary_degree'],row['depth'],row['term'],
1324 row['partial'],row['error'],'infinity' if row['valuation'] is None else row['valuation']])
1325 with paths[4].open('w',encoding='utf-8',newline='') as stream:
1326 writer=csv.writer(stream,lineterminator='\n')
1327 writer.writerow(['transform','chart','row','column','value','coordinates'])
1328 def write_matrix(name,chart,matrix,coordinates):
1329 for n,row in enumerate(matrix):
1330 for j,value in enumerate(row):
1331 writer.writerow([name,chart,n,j,str(value),coordinates])
1332 for chart in (2,3,5,7):
1333 for name,kernel in [('local_forward',local_kernel(chart,12)),
1334 ('local_inverse',local_inverse_kernel(chart,12))]:
1335 write_matrix(name,chart,coefficient_matrix(kernel,12),'even EGF numerators')
1336 write_matrix('pascal_reconstruction','',pascal_matrix(12),'sequence from differences')
1337 write_matrix('pascal_difference','',pascal_matrix(12,inverse=True),'differences from sequence')
1338 write_matrix('sampling_M2_s0','',report['matrix_examples']['sampling_two'],'Newton differences')
1339 write_matrix('source_base_change','3_to_5',report['base_change_examples']['three_to_five'],'even EGF source numerators')
1340 write_matrix('source_base_change','3_to_2',report['base_change_examples']['three_to_two'],'even EGF source numerators')
1341 with paths[5].open('w',encoding='utf-8',newline='') as stream:
1342 writer=csv.writer(stream,lineterminator='\n')
1343 writer.writerow(['chart','depth','proved_modulus','partial','residue','error','exact_error_valuation'])
1344 for row in report['base_change_examples']['congruence_ladders']:
1345 writer.writerow([row[key] for key in ('chart','depth','modulus','partial','residue','error','valuation')])
1346 with paths[6].open('w',encoding='utf-8',newline='') as stream:
1347 writer=csv.writer(stream,lineterminator='\n')
1348 writer.writerow(['transform','row','column','value'])
1349 for name in ('polynomial_connection','polynomial_inverse','recurrence_analysis',
1350 'recurrence_synthesis','sequence_transfer','sequence_inverse'):
1351 for i,row in enumerate(report['fibonacci'][name]):
1352 for j,value in enumerate(row):
1353 writer.writerow([name,i,j,value])
1354 with paths[7].open('w',encoding='utf-8',newline='') as stream:
1355 writer=csv.writer(stream,lineterminator='\n')
1356 keys=['step','refinement','depth','target','partial','error','modulus','valuation']
1357 writer.writerow(keys)
1358 for row in report['literature']['precision']:
1359 writer.writerow(['infinity' if row[key] is None else row[key] for key in keys])
1360 with paths[8].open('w',encoding='utf-8',newline='') as stream:
1361 writer=csv.writer(stream,lineterminator='\n')
1362 writer.writerow(['source_degree','operator','row','column','value'])
1363 for degree,record in report['literature']['calculus'].items():
1364 for name in ('integral','derivative','euler_derivative','connection','inverse_connection'):
1365 for i,row in enumerate(record[name]):
1366 for j,value in enumerate(row):
1367 writer.writerow([degree,name,i,j,value])
1368 with paths[9].open('w',encoding='utf-8',newline='') as stream:
1369 writer=csv.writer(stream,lineterminator='\n')
1370 writer.writerow(['address','index','normalized_C','original_U','newton_coefficient','modulus','modular_C','phi_mod_256'])
1371 atlas=report['atlas']
1372 for address in range(4):
1373 for n in range(256):
1374 exact=atlas['quartic'][address][n] if n<48 else ''
1375 original=str(factorial(address)*F(exact)) if exact else ''
1376 newton=atlas['newton_prefix'][address][n] if n<13 else ''
1377 writer.writerow([address,n,exact,original,newton,8192,atlas['modular_coefficients'][address][n],atlas['phi_residues'][address][n]])
1378 with paths[10].open('w',encoding='utf-8',newline='') as stream:
1379 writer=csv.DictWriter(stream,fieldnames=['bits','modulus','low','low_image','high','high_image','target','chosen'],lineterminator='\n')
1380 writer.writeheader()
1381 writer.writerows(report['atlas']['inverse_trace'])
1382 with paths[11].open('w',encoding='utf-8',newline='') as stream:
1383 writer=csv.writer(stream,lineterminator='\n')
1384 writer.writerow(['full_degree','signed_euler','logarithmic_hurwitz','logarithmic_ordinary','mod_3','mod_5','mod_7'])
1385 for n in range(25):
1386 writer.writerow([n,report['atlas']['full_euler'][n],
1387 report['atlas']['logarithmic_euler'][n] if n<17 else '',
1388 report['atlas']['logarithmic_ordinary'][n] if n<17 else '',
1389 *[row['values'][n] for row in report['atlas']['modes']]])
1390 with paths[12].open('w',encoding='utf-8',newline='') as stream:
1391 keys=['r','N','x','retained_depth','partial','tail','H','H_valuation','tail_valuation','binomial_valuation','minimum_summand_valuation']
1392 writer=csv.writer(stream,lineterminator='\n')
1393 writer.writerow(keys)
1394 for row in report['edge']['rows']:writer.writerow(['infinity' if row[k] is None else row[k] for k in keys])
1395
1396
1397 def main():
1398 parser = argparse.ArgumentParser(description=__doc__)
1399 parser.add_argument('--output',type=Path,help='Save JSON and CSV results in this directory')
1400 args = parser.parse_args()
1401 report = compute()
1402 if args.output:
1403 save_results(args.output,report)
1404 print('Secant numbers:', ', '.join(map(str,report['secant'][:7])))
1405 print('Quartic reciprocal:', ', '.join(map(str,report['quartic'][:5])))
1406 print('Odd-length and dyadic reconstructions: E_8 = 1385')
1407 print('Dyadic E_8 tail valuations: 5, 8, 13, 15, exact')
1408 print('Matrix checks: inverse, composition, factorial normalization, packets, Pascal sampling')
1409 print('Sampling cutoff control: 3 columns give 32; 5 columns give 1376')
1410 print('Base changes: exact cocycle; prime-dependent precision checked')
1411 print('Supercongruence combination: 8 mod 81 and 9 mod 32 give 1385 mod 2592')
1412 print('Fibonacci polynomial coordinates: E_4(x) = 2 B_0 + 5 B_1 - 3 B_2 - 2 B_3 + B_4')
1413 print('Fibonacci recurrence defects: 5, -65, 1441, -51845')
1414 print('Reversible sequence transfer: 1, -1, 5, -61, 1385 -> 1, 1, 2, 3, 5')
1415 print('Paper IX Fibonacci packets: step 3 gives u[k+2] = 4*u[k+1] + u[k]')
1416 print('Fibonacci precision ladder: F_37 = 24157817; tail valuations 3, 6, 9, exact')
1417 print('Euler/Fibonacci calculus: full round trip is I_4; integral = 5/4')
1418 print('Quartic address gains: normalized 1, 1, 4, 2; original 1, 1, 5, 3')
1419 print('Dyadic inverse: address 2, target 13 mod 64 -> index',report['atlas']['inverse_index'])
1420 print('Finite mode inverse mod 5: weights 1, 1, -1, -1, 1')
1421 print('Logarithmic Euler coordinates: 1, 0, -1, 3, -6, 0, 90, -630, 2520')
1422 print('Dyadic edge: exact tail factorization and sparse-binary row valuations checked')
1423 print('Signed Euler difference coordinates:', ', '.join(map(str,report['difference_coordinates'])))
1424 print('Spacing-two lane: 1, 5, 1385; differences: 1, 4, 1376')
1425 print('Refinement: 2*4 + 1376 = 1384')
1426 print('Prime residual: 2*(-2) + 8 = 4')
1427 print('Newton residual orders:', ', '.join(map(str,report['newton_residual_orders'])))
1428 print('PASS: Euler packets, polynomials and residuals')
1429
1430
1431 if __name__ == '__main__':
1432 main()