Roots of unity, projectors, and falling factorials / exact Python
Reconstruct x^5 from three transported character packets. Includes exact Newton/Stirling inverses, difference-transfer identities, a six-packet refinement, an induced star product, direct function enumeration, the spectral number operator, filtered EGFs, translation carry blocks, finite-characteristic obstructions, and counterexamples. No complex floating-point arithmetic or third-party package is required.
Read the complete mathematical derivation · Download standalone Python · Source and results ZIP · Source identity record
Exact generated files: newton-stirling.csv · image-histograms.csv · packet-egfs.csv · spectral-projectors.csv · translation-blocks.csv · shifted-packet-ledger.csv · finite-characteristic-obstruction.csv · packet-contributions.csv · difference-transfer.csv · refinement.csv · projector-matrix.csv · projection-diagram.svg · function-image-histogram.svg · translation-carry.svg · results.json
Run python -B roots_unity_falling_factorial.py to check, or use --output roots-results for the CSV and SVG outputs.
1 #!/usr/bin/env python3
2 """Exact roots-of-unity, Newton and falling-factorial projector example.
3
4 All computation uses rational arithmetic in the Python standard library.
5 Roots of unity are handled through the equivalent exact cyclic trace,
6 not floating-point evaluations. These are classical finite identities.
7 """
8 from __future__ import annotations
9
10 import argparse
11 import csv
12 import json
13 from fractions import Fraction
14 from itertools import product
15 from math import comb, factorial
16 from pathlib import Path
17
18 Q = Fraction
19 MAX_DEGREE = 8
20 ROW_FIVE = (0, 1, 15, 25, 10, 1)
21 FILES = (
22 "newton-stirling.csv",
23 "image-histograms.csv",
24 "packet-egfs.csv",
25 "spectral-projectors.csv",
26 "translation-blocks.csv",
27 "shifted-packet-ledger.csv",
28 "finite-characteristic-obstruction.csv",
29 "packet-contributions.csv",
30 "difference-transfer.csv",
31 "refinement.csv",
32 "projector-matrix.csv",
33 "projection-diagram.svg",
34 "function-image-histogram.svg",
35 "translation-carry.svg",
36 "results.json",
37 )
38
39
40 def trim(values):
41 """Canonical ascending-power coefficient tuple."""
42 out = [Q(x) for x in values]
43 while len(out) > 1 and out[-1] == 0:
44 out.pop()
45 return tuple(out) if out else (Q(0),)
46
47
48 def add(a, b):
49 out = [Q(0)] * max(len(a), len(b))
50 for i, x in enumerate(a):
51 out[i] += x
52 for i, x in enumerate(b):
53 out[i] += x
54 return trim(out)
55
56
57 def scale(a, c):
58 return trim(Q(c) * v for v in a)
59
60
61 def multiply(a, b):
62 out = [Q(0)] * (len(a) + len(b) - 1)
63 for i, x in enumerate(a):
64 for j, y in enumerate(b):
65 out[i + j] += x * y
66 return trim(out)
67
68
69 def evaluate(a, x):
70 value = Q(0)
71 for coefficient in reversed(a):
72 value = value * x + coefficient
73 return value
74
75
76 def power(n):
77 if n < 0:
78 raise ValueError("Degree must be nonnegative")
79 return (Q(0),) * n + (Q(1),)
80
81
82 def falling(k):
83 """x_(k)=x(x-1)...(x-k+1), with x_(0)=1."""
84 if k < 0:
85 raise ValueError("Falling-factorial degree must be nonnegative")
86 result = (Q(1),)
87 for j in range(k):
88 result = multiply(result, (Q(-j), Q(1)))
89 return result
90
91
92 def shift(a):
93 """Power-basis coefficients of f(x+1)."""
94 out = [Q(0)] * len(a)
95 for k, coefficient in enumerate(a):
96 for j in range(k + 1):
97 out[j] += coefficient * comb(k, j)
98 return trim(out)
99
100
101 def delta(a):
102 """Forward difference, f(x+1)-f(x)."""
103 return add(shift(a), scale(a, -1))
104
105
106 def stirling_second(n):
107 rows = [(1,)]
108 for i in range(1, n + 1):
109 prev = rows[-1]
110 row = [0] * (i + 1)
111 for k in range(1, i + 1):
112 row[k] = (prev[k-1] if k-1 < len(prev) else 0) + (
113 k * prev[k] if k < len(prev) else 0)
114 rows.append(tuple(row))
115 return tuple(rows)
116
117
118 def signed_stirling_first(n):
119 rows = [(1,)]
120 for i in range(1, n + 1):
121 prev = rows[-1]
122 row = [0] * (i + 1)
123 for k in range(i + 1):
124 row[k] = ((prev[k-1] if k else 0) -
125 (i - 1) * (prev[k] if k < len(prev) else 0))
126 rows.append(tuple(row))
127 return tuple(rows)
128
129
130 def from_newton(coefficients):
131 """Phi(sum a_k z^k)=sum a_k x_(k). Linear, not multiplicative."""
132 out = (Q(0),)
133 for k, coefficient in enumerate(coefficients):
134 out = add(out, scale(falling(k), coefficient))
135 return out
136
137
138 def to_newton(poly):
139 """Phi inverse via Newton interpolation at 0,1,...,degree."""
140 poly = trim(poly)
141 degree = len(poly) - 1
142 diffs = [evaluate(poly, j) for j in range(degree + 1)]
143 result = []
144 for k in range(degree + 1):
145 result.append(diffs[0] / factorial(k))
146 diffs = [diffs[j+1] - diffs[j] for j in range(len(diffs)-1)]
147 return tuple(result)
148
149
150 def project_newton(poly, modulus, residue):
151 """Transport a character projector through Phi to polynomial x."""
152 if not isinstance(modulus, int) or modulus < 2:
153 raise ValueError("Require integer modulus at least two")
154 if not isinstance(residue, int) or not 0 <= residue < modulus:
155 raise ValueError("Residue must be in 0,...,modulus-1")
156 coordinates = to_newton(poly)
157 return from_newton(tuple(c if k % modulus == residue else Q(0)
158 for k, c in enumerate(coordinates)))
159
160
161 def naive_power_project(poly, modulus, residue):
162 """Valid ordinary-power projector; not the transported Newton projector."""
163 if type(modulus) is not int or modulus < 2:
164 raise ValueError("Require integer modulus at least two")
165 if type(residue) is not int or not 0 <= residue < modulus:
166 raise ValueError("Residue must be in 0,...,modulus-1")
167 return trim(v if k % modulus == residue else Q(0)
168 for k, v in enumerate(poly))
169
170
171 def trace_character_selector(modulus, k, residue):
172 """(1/m)Tr(C^(k-r)) for the regular cyclic permutation representation.
173
174 C has the m eigenvalues of m-th roots of unity; the trace is exactly
175 m or 0. This checks Fourier orthogonality using integers only.
176 """
177 if modulus < 2 or not 0 <= residue < modulus:
178 raise ValueError("Invalid cyclic character filter")
179 displacement = k - residue
180 fixed_points = sum(i == (i + displacement) % modulus
181 for i in range(modulus))
182 return Q(fixed_points, modulus)
183
184
185 def star(a, b):
186 """Transport ordinary multiplication of Newton-coordinate polynomials."""
187 ac, bc = to_newton(a), to_newton(b)
188 result = [Q(0)] * (len(ac) + len(bc) - 1)
189 for i, x in enumerate(ac):
190 for j, y in enumerate(bc):
191 result[i+j] += x * y
192 return from_newton(result)
193
194
195 def shift_newton_derivative(coordinates):
196 """Derivative on z coordinates, corresponding to Delta on x."""
197 return tuple(Q(k+1)*coordinates[k+1]
198 for k in range(len(coordinates)-1)) or (Q(0),)
199
200
201
202 def translate(a, offset):
203 """f(x+offset), computed independently by binomial expansion."""
204 out = [Q(0)] * len(a)
205 for k, coefficient in enumerate(a):
206 for j in range(k + 1):
207 out[j] += Q(coefficient) * comb(k, j) * Q(offset)**(k-j)
208 return trim(out)
209
210
211
212 def translate_via_newton(poly, offset):
213 """Independent translation from the falling-factorial binomial identity.
214
215 (x+h)_(k) = sum_{j<=k} choose(k,j)*(h)_(k-j)*(x)_(j).
216 Valid for any rational shift h; not based on power-basis substitution.
217 """
218 h = Q(offset)
219 coefficients = to_newton(poly)
220 out = [Q(0)] * len(coefficients)
221 for k, value in enumerate(coefficients):
222 for j in range(k+1):
223 out[j] += value * comb(k,j) * evaluate(falling(k-j),h)
224 return from_newton(out)
225
226
227 def translate_packet_block(poly, modulus, source, target, offset):
228 """Q_target T_h Q_source in falling-factorial coordinates.
229
230 Constructed from the explicit binomial-type connection coefficients,
231 independently of translate() and project_newton().
232 """
233 if type(modulus) is not int or modulus < 2:
234 raise ValueError("Require integer modulus >= 2")
235 if type(source) is not int or not 0 <= source < modulus:
236 raise ValueError("Invalid source residue")
237 if type(target) is not int or not 0 <= target < modulus:
238 raise ValueError("Invalid target residue")
239 h = Q(offset)
240 coefficients = to_newton(poly)
241 out = [Q(0)] * len(coefficients)
242 for k,value in enumerate(coefficients):
243 if k % modulus != source:
244 continue
245 for j in range(k+1):
246 if j % modulus == target:
247 out[j] += value * comb(k,j) * evaluate(falling(k-j),h)
248 return from_newton(out)
249
250
251 def mod_polynomial(poly, modulus):
252 """Integer-coefficient polynomial reduced modulo an integer >= 2.
253
254 This is deliberately a bounded exact finite-characteristic witness, not
255 an implementation of general division or spectral functional calculus.
256 """
257 if type(modulus) is not int or modulus < 2:
258 raise ValueError("Require integer modulus >= 2")
259 entries = tuple(Q(v) for v in poly)
260 if any(v.denominator != 1 for v in entries):
261 raise ValueError("Require integral polynomial coefficients")
262 values = [v.numerator % modulus for v in entries]
263 while len(values)>1 and values[-1]==0:
264 values.pop()
265 return tuple(values) if values else (0,)
266
267
268 def backward_difference(a):
269 """Nabla f(x) = f(x) - f(x-1)."""
270 return add(a, scale(translate(a, -1), -1))
271
272
273 def newton_degree_operator(poly):
274 """Number operator N = x*Nabla: N(x_(k))=k*x_(k)."""
275 return multiply(power(1), backward_difference(poly))
276
277
278 def spectral_project_newton(poly, modulus, residue):
279 """Exact projector as rational Lagrange polynomial in N on Q[x]_{<=d}."""
280 if not isinstance(modulus, int) or modulus < 2:
281 raise ValueError("Require integer modulus at least two")
282 if not isinstance(residue, int) or not 0 <= residue < modulus:
283 raise ValueError("Residue must be in 0,...,modulus-1")
284 poly = trim(poly)
285 degree = len(poly) - 1
286 result = (Q(0),)
287 for k in range(degree + 1):
288 if k % modulus != residue:
289 continue
290 candidate = poly
291 for j in range(degree + 1):
292 if j != k:
293 candidate = scale(
294 add(newton_degree_operator(candidate),scale(candidate, -j)),
295 Q(1, k-j))
296 result = add(result, candidate)
297 return result
298
299
300 def image_count_histogram(n, codomain):
301 """Independent exhaustive enumeration of functions [n] -> [codomain].
302
303 For n=0, the empty function is counted once, even for empty codomain.
304 No Stirling recurrence, Newton transform or roots of unity is consulted.
305 """
306 if type(n) is not int or type(codomain) is not int or n < 0 or codomain < 0:
307 raise ValueError("n and codomain must be nonnegative integers")
308 counts = [0] * (n+1)
309 for image in product(range(codomain), repeat=n):
310 counts[len(set(image))] += 1
311 return tuple(counts)
312
313
314 def series_product_truncated(a, b, degree):
315 out = [Q(0)] * (degree + 1)
316 for i, value in enumerate(a):
317 for j, other in enumerate(b):
318 if i + j <= degree:
319 out[i+j] += value * other
320 return tuple(out)
321
322
323 def packet_egf_coefficients(x, modulus, residue, degree):
324 """[t^n] of sum_{k congruent residue} binom(x,k)*(exp(t)-1)^k.
325
326 Formal power series, with rational x; no evaluation at complex roots.
327 """
328 if not isinstance(modulus, int) or modulus < 2:
329 raise ValueError("Require integer modulus at least two")
330 if not isinstance(residue, int) or not 0 <= residue < modulus:
331 raise ValueError("Residue must be in 0,...,modulus-1")
332 if type(degree) is not int or degree < 0:
333 raise ValueError("Degree must be a nonnegative integer")
334 x = Q(x)
335 u = (Q(0),) + tuple(Q(1, factorial(k)) for k in range(1, degree+1))
336 power_u = (Q(1),) + (Q(0),)*degree
337 result = [Q(0)]*(degree+1)
338 for k in range(degree+1):
339 if k % modulus == residue:
340 weight = evaluate(falling(k), x) / factorial(k)
341 for n in range(degree+1):
342 result[n] += weight*power_u[n]
343 power_u = series_product_truncated(power_u,u,degree)
344 return tuple(result)
345
346
347 def example_row():
348 row = stirling_second(5)[5]
349 polynomial = power(5)
350 pieces = tuple(project_newton(polynomial, 3, r) for r in range(3))
351 return row, polynomial, pieces
352
353
354 def check():
355 counts = {}
356
357 def verify(cond, category, reason):
358 if not cond:
359 raise AssertionError(category + ": " + reason)
360 counts[category] = counts.get(category, 0) + 1
361
362 row, p, packets = example_row()
363 verify(row == ROW_FIVE, "anchor", "Stirling fifth row")
364 verify(to_newton(p) == tuple(map(Q, ROW_FIVE)), "anchor", "Newton fifth row")
365 verify(from_newton(ROW_FIVE) == p, "anchor", "recover x^5")
366 verify(tuple(evaluate(a, 4) for a in packets) == (600, 244, 180),
367 "anchor", "three exact packet values")
368 verify(sum(evaluate(a, 4) for a in packets) == 1024,
369 "anchor", "600+244+180=4^5")
370 verify(tuple(evaluate(delta(a), 4) for a in packets) == (900, 961, 240),
371 "anchor", "difference packet values")
372 verify(sum(evaluate(delta(a), 4) for a in packets) == 2101,
373 "anchor", "5^5-4^5")
374 verify(evaluate(naive_power_project(p, 3, 2), 4) == 1024,
375 "anchor", "incorrect raw-degree filter")
376 verify(evaluate(packets[2], 4) == 180, "anchor", "correct transferred filter")
377 verify(star(power(1), power(1)) == falling(2),
378 "anchor", "star multiplication on first-order coordinates")
379 verify(multiply(power(1), power(1)) != star(power(1), power(1)),
380 "anchor", "Phi does not preserve ordinary products")
381 verify(delta(power(2)) == (Q(1), Q(2)),
382 "anchor", "finite-difference Leibniz obstruction")
383
384 second, first = stirling_second(MAX_DEGREE), signed_stirling_first(MAX_DEGREE)
385 for n in range(MAX_DEGREE + 1):
386 recovered = to_newton(power(n))
387 for k in range(n+1):
388 verify(recovered[k] == second[n][k], "stirling",
389 f"independent difference triangle n={n},k={k}")
390 verify(falling(k)[0] == int(k == 0), "stirling",
391 f"falling zero coefficient k={k}")
392 for k in range(n+1):
393 verify(sum(first[n][j] * second[j][k] for j in range(k, n+1))
394 == int(n == k), "stirling", f"inverse first*second {n},{k}")
395 verify(sum(second[n][j] * first[j][k] for j in range(k, n+1))
396 == int(n == k), "stirling", f"inverse second*first {n},{k}")
397 verify(from_newton(recovered) == power(n), "stirling",
398 f"Newton roundtrip x^{n}")
399 verify(falling(n) == tuple(map(Q, first[n])), "stirling",
400 f"signed first-kind polynomial n={n}")
401
402 sample_polys = [power(i) for i in range(MAX_DEGREE+1)]
403 sample_polys += [(Q(2),Q(-3),Q(0),Q(5)), (Q(1,2),Q(0),Q(-4),Q(1),Q(2))]
404 for poly in sample_polys:
405 verify(from_newton(to_newton(poly)) == trim(poly),
406 "roundtrip", "polynomial Newton roundtrip")
407 for modulus in range(2,7):
408 pieces = [project_newton(poly, modulus, r) for r in range(modulus)]
409 whole = (Q(0),)
410 for part in pieces:
411 whole = add(whole, part)
412 verify(whole == trim(poly), "projection", "complete residue family")
413 for r in range(modulus):
414 verify(project_newton(pieces[r], modulus, r) == pieces[r],
415 "projection", "idempotent")
416 verify(project_newton(delta(pieces[r]), modulus,(r-1)%modulus)
417 == delta(pieces[r]), "difference",
418 "Delta Q_r = Q_(r-1) Delta")
419 verify(delta(pieces[r]) ==
420 project_newton(delta(poly), modulus,(r-1)%modulus),
421 "difference", "strong differential transfer")
422 for j in range(modulus):
423 if r != j:
424 verify(project_newton(pieces[r], modulus, j) == (Q(0),),
425 "projection", "pairwise annihilation")
426 verify(trace_character_selector(modulus,r,r)==1,
427 "fourier", "character self-sum")
428 verify(delta(poly) ==
429 from_newton(shift_newton_derivative(to_newton(poly))),
430 "difference", "Delta = Phi D Phi inverse")
431 for k in range(len(poly)):
432 for r in range(modulus):
433 verify(trace_character_selector(modulus,k,r) ==
434 int(k % modulus == r), "fourier",
435 "cyclic trace equals roots-of-unity average")
436 for factor in (2,3):
437 refined=modulus*factor
438 for r in range(modulus):
439 children=(Q(0),)
440 for q in range(factor):
441 children=add(children,project_newton(poly,refined,r+modulus*q))
442 verify(children==pieces[r],"refinement",
443 "parent equals disjoint residue children")
444 # Independent power-basis projectors are legitimate, but are not the
445 # Newton-transferred projectors unless the chosen coordinate change
446 # preserves those specific degree labels.
447 for poly in sample_polys:
448 for modulus in range(2,7):
449 power_parts=[naive_power_project(poly,modulus,r)
450 for r in range(modulus)]
451 reconstructed=(Q(0),)
452 for part in power_parts:
453 reconstructed=add(reconstructed,part)
454 verify(reconstructed==trim(poly),"power_projector",
455 "ordinary power-degree filters reconstruct the polynomial")
456 for r in range(modulus):
457 for target in range(modulus):
458 verify(naive_power_project(power_parts[r],modulus,target)
459 ==(power_parts[r] if r==target else (Q(0),)),
460 "power_projector",
461 "ordinary power-degree filters are orthogonal projectors")
462
463 # A second, *independent* interpretation: count actual functions by
464 # the number of different labels in their image, without using S(n,k).
465 verify(image_count_histogram(5,4)==(0,4,180,600,240,0),
466 "enumeration","all 1,024 functions [5] -> [4] by image size")
467 for codomain in range(5):
468 for n in range(7):
469 observed=image_count_histogram(n,codomain)
470 verify(sum(observed)==codomain**n,"enumeration",
471 "all functions counted, including the empty function")
472 for k in range(n+1):
473 expected=second[n][k]*evaluate(falling(k),codomain)
474 verify(observed[k]==expected,"enumeration",
475 "independent image enumeration vs falling Stirling count")
476
477 # The finite-difference number operator N=x*(I-E^{-1}) has eigenvectors
478 # (x)_(k) with eigenvalues k. A finite degree cutoff makes its spectral
479 # projectors rational Lagrange polynomials in N, no roots needed.
480 for k in range(MAX_DEGREE+1):
481 verify(newton_degree_operator(falling(k))==scale(falling(k),k),
482 "spectral","falling factorial eigenvalue")
483 verify(newton_degree_operator(power(2))==(Q(0),Q(-1),Q(2)),
484 "spectral","N x^2=2x^2-x, not ordinary Euler derivative")
485 for polynomial in sample_polys:
486 verify(add(newton_degree_operator(delta(polynomial)),
487 scale(delta(newton_degree_operator(polynomial)),-1))
488 ==scale(delta(polynomial),-1),"spectral",
489 "[N,Delta] = -Delta")
490 for n in range(7):
491 for modulus in (2,3,4,5):
492 for r in range(modulus):
493 verify(spectral_project_newton(power(n),modulus,r)
494 ==project_newton(power(n),modulus,r),
495 "spectral","spectral Lagrange projector equals conjugation")
496
497 # A third independent route via the EGF of functions with image size
498 # in a specified residue class. Compare all formal t coefficients.
499 for x in (Q(0),Q(1),Q(4),Q(1,2),Q(-3,2)):
500 for modulus in (2,3,4):
501 egfs=[packet_egf_coefficients(x,modulus,r,MAX_DEGREE)
502 for r in range(modulus)]
503 for n in range(MAX_DEGREE+1):
504 verify(sum(egfs[r][n] for r in range(modulus)) ==
505 Q(x**n,factorial(n)),"packet_egf",
506 "sum of filtered EGFs is exp(x*t)")
507 for r in range(modulus):
508 verify(egfs[r][n] ==
509 Q(evaluate(project_newton(power(n),modulus,r),x),
510 factorial(n)),"packet_egf",
511 "EGF coefficients reproduce transported projection")
512 verify(tuple(packet_egf_coefficients(4,3,r,5)[5]*factorial(5)
513 for r in range(3))==(600,244,180),
514 "packet_egf","fifth EGF coefficient matches image partition")
515
516 # Translation, filtered packets and all cross-channel carry blocks.
517 source_packets=tuple(project_newton(p,3,r) for r in range(3))
518 shifted_target=tuple(project_newton(translate(p,1),3,r) for r in range(3))
519 verify(tuple(evaluate(q,4) for q in shifted_target)==(1561,484,1080),
520 "translation_anchor","shifted output packet values")
521 verify(tuple(evaluate(translate(q,1),4) for q in source_packets)
522 ==(1500,1205,420),"translation_anchor",
523 "individually shifted source packets")
524 verify(tuple(evaluate(delta(q),4) for q in source_packets)==(900,961,240),
525 "translation_anchor","individual carry terms")
526 for r in range(3):
527 verify(shifted_target[r] ==
528 add(source_packets[r],delta(source_packets[(r+1)%3])),
529 "translation_anchor","Q_r T_1 = Q_r + Delta Q_(r+1)")
530 verify(sum(evaluate(q,4) for q in shifted_target)==3125,
531 "translation_anchor","shifted packets recover fifth power")
532 verify(shifted_target!=tuple(translate(q,1) for q in source_packets),
533 "translation_anchor","translation does not preserve original packet labels")
534
535 for poly in sample_polys:
536 for h in (Q(-1),Q(0),Q(1),Q(2),Q(3),Q(1,2)):
537 shifted=translate(poly,h)
538 verify(shifted==translate_via_newton(poly,h),"translation",
539 "Newton binomial translation equals independent power translation")
540 for modulus in range(2,7):
541 source_parts=[project_newton(poly,modulus,r)
542 for r in range(modulus)]
543 target_parts=[project_newton(shifted,modulus,r)
544 for r in range(modulus)]
545 for r in range(modulus):
546 aggregate=(Q(0),)
547 for source in range(modulus):
548 block=translate_packet_block(poly,modulus,source,r,h)
549 verify(block == project_newton(
550 translate(source_parts[source],h),modulus,r),
551 "translation_blocks","independent cross-packet transfer")
552 aggregate=add(aggregate,block)
553 if h.denominator == 1 and h.numerator>=0:
554 lag=(source-r)%modulus
555 # For nonnegative integer shift h, the falling
556 # binomial identity has no terms of lag > h.
557 if all((k-j)>h for k in range(len(poly))
558 for j in range(k+1)
559 if k%modulus==source and j%modulus==r):
560 verify(block==(Q(0),),"translation_blocks",
561 "zero outside integer-shift bandwidth")
562 verify(target_parts[r]==aggregate,"translation",
563 "all packet blocks reconstruct target packet")
564 if h==1:
565 for source in range(modulus):
566 for r in range(modulus):
567 block=translate_packet_block(poly,modulus,source,r,h)
568 rhs=(Q(0),)
569 if r==source:
570 rhs=add(rhs,source_parts[source])
571 if r==(source-1)%modulus:
572 rhs=add(rhs,delta(source_parts[source]))
573 verify(block==rhs,"translation_blocks",
574 "one-step transfer is cyclic two-tap")
575 if h==0:
576 for source in range(modulus):
577 for r in range(modulus):
578 verify(translate_packet_block(
579 poly,modulus,source,r,h)==
580 (source_parts[source] if r==source else (Q(0),)),
581 "translation_blocks","zero shift gives diagonal transfer")
582
583 # The signed first/second Stirling matrices are integral and mutually
584 # inverse. Unlike Newton's divided-difference formula, the change of
585 # basis remains invertible after reduction modulo any prime.
586 for prime in (2,3,5):
587 for n in range(MAX_DEGREE+1):
588 for k in range(n+1):
589 row_first_second=sum(first[n][j]*second[j][k]
590 for j in range(k,n+1))
591 row_second_first=sum(second[n][j]*first[j][k]
592 for j in range(k,n+1))
593 verify(row_first_second%prime == int(n==k),
594 "modular_basis","signed-first times second inverse mod p")
595 verify(row_second_first%prime == int(n==k),
596 "modular_basis","second times signed-first inverse mod p")
597
598 # Over F3 the cubic forward difference vanishes on every polynomial;
599 # x^3-x is nevertheless a nonzero polynomial despite being zero as
600 # a function on the three-point field.
601 for n in range(MAX_DEGREE+1):
602 poly=power(n)
603 cubed=delta(delta(delta(poly)))
604 verify(mod_polynomial(cubed,3)==(0,),
605 "modular_basis","third forward difference vanishes modulo 3")
606 verify(mod_polynomial(translate(poly,3),3)==mod_polynomial(poly,3),
607 "modular_basis","shift by characteristic is polynomial identity")
608 verify(mod_polynomial(falling(3),3)!=(0,),"modular_basis",
609 "x^3-x is nonzero as a polynomial over F3")
610 verify(all(evaluate(falling(3),a).numerator%3==0 for a in range(3)),
611 "modular_basis","same polynomial is zero as a function on F3")
612
613 # The same carry in the function-image generating functions. This checks
614 # the binomial coefficient identity in an independent formal-series route.
615 u=(Q(0),)+tuple(Q(1,factorial(j)) for j in range(1,MAX_DEGREE+1))
616 for x in (Q(0),Q(1),Q(4),Q(1,2),Q(-3,2)):
617 for modulus in (2,3,4):
618 for residue in range(modulus):
619 before=packet_egf_coefficients(x,modulus,residue,MAX_DEGREE)
620 after=packet_egf_coefficients(x+1,modulus,residue,MAX_DEGREE)
621 prior=packet_egf_coefficients(x,modulus,(residue-1)%modulus,
622 MAX_DEGREE)
623 carry=series_product_truncated(u,prior,MAX_DEGREE)
624 for n in range(MAX_DEGREE+1):
625 verify(after[n]==before[n]+carry[n],"egf_carry",
626 "G_r(t,x+1)=G_r(t,x)+(exp(t)-1)G_(r-1)(t,x)")
627
628 # Characteristic-three obstruction to an unjustified spectral lift:
629 # falling factorials remain an integral unitriangular basis; however
630 # number-operator eigenvalues 0 and 3 collide mod 3.
631 one=falling(0)
632 triple=falling(3)
633 verify(mod_polynomial(triple,3)==(0,2,0,1),"characteristic",
634 "falling degree three is a nonzero polynomial in F3[x]")
635 verify(mod_polynomial(newton_degree_operator(one),3)==(0,),
636 "characteristic","eigenvalue zero on the constant polynomial")
637 verify(mod_polynomial(newton_degree_operator(triple),3)==(0,),
638 "characteristic","eigenvalue three collapses to zero in F3")
639 verify(mod_polynomial(project_newton(one,2,0),3)==(1,),
640 "characteristic","even character projector keeps degree zero")
641 verify(mod_polynomial(project_newton(triple,2,0),3)==(0,),
642 "characteristic","even character projector removes degree three")
643 verify(mod_polynomial(project_newton(triple,2,1),3)==(0,2,0,1),
644 "characteristic","odd character projector keeps degree three")
645 verify(mod_polynomial(falling(3),3)!=mod_polynomial(falling(0),3),
646 "characteristic","eigenvectors remain distinct")
647 verify(mod_polynomial(translate(triple,3),3)==mod_polynomial(triple,3),
648 "characteristic","shift by characteristic is identity on F3[x]")
649
650 # All entries of Q in the ordinary-power basis, computed independently
651 # from signed Stirling and the second-kind triangle.
652 for degree in range(0,7):
653 for modulus in (2,3,4,5):
654 for r in range(modulus):
655 for j in range(degree+1):
656 projected=project_newton(power(j),modulus,r)
657 for i in range(degree+1):
658 actual=projected[i] if i<len(projected) else Q(0)
659 expected=sum(first[k][i]*second[j][k]
660 for k in range(i,j+1) if k%modulus==r)
661 verify(actual==expected,"matrix",
662 "Q=M diag_r M^-1 in monomial coordinates")
663
664 examples=[
665 (power(1),power(1)),
666 ((Q(1),Q(2),Q(1)),(Q(3),Q(-1),Q(0),Q(1))),
667 ((Q(-2),Q(1),Q(0),Q(3)),(Q(2),Q(1),Q(-1))),
668 ]
669 for f,g in examples:
670 verify(to_newton(star(f,g)) ==
671 trim(tuple(sum(
672 to_newton(f)[i]*to_newton(g)[k-i]
673 for i in range(len(to_newton(f)))
674 if 0<=k-i<len(to_newton(g)))
675 for k in range(len(to_newton(f))+len(to_newton(g))-1))),
676 "star","Newton multiplication law")
677 verify(delta(star(f,g)) ==
678 add(star(delta(f),g),star(f,delta(g))),
679 "star","Delta is star derivation")
680 # The ordinary discrete Leibniz rule has a shifted first factor.
681 verify(delta(multiply(f,g)) ==
682 add(multiply(delta(f),g),
683 multiply(shift(f),delta(g))),
684 "star","correct shifted ordinary Leibniz rule")
685 for m in (2,3,4):
686 for r in range(m):
687 right=(Q(0),)
688 for a in range(m):
689 right=add(right,star(project_newton(f,m,a),
690 project_newton(g,m,(r-a)%m)))
691 verify(project_newton(star(f,g),m,r)==right,
692 "star","cyclic packet multiplication")
693 for bad in (0,1,-2,2.5):
694 try:
695 project_newton(p,bad,0)
696 except ValueError:
697 verify(True,"boundaries","reject invalid modulus")
698 else:
699 verify(False,"boundaries","invalid modulus accepted")
700 for bad in (-1,3,1.5):
701 try:
702 project_newton(p,3,bad)
703 except ValueError:
704 verify(True,"boundaries","reject invalid residue")
705 else:
706 verify(False,"boundaries","invalid residue accepted")
707 return counts
708
709
710 def svg_diagram():
711 """Accessible, vector-only illustration generated from the same exact data."""
712 row, poly, packets = example_row()
713 values = [int(evaluate(p,4)) for p in packets]
714 pieces = ["25(x)_3", "(x)_1 + 10(x)_4", "15(x)_2 + (x)_5"]
715 fills = ["#d4f1e8","#e0e9ff","#fff0d9"]
716 strokes = ["#156b5b","#365d9d","#9c5816"]
717 out=[
718 '<svg xmlns="http://www.w3.org/2000/svg" width="1080" height="580" viewBox="0 0 1080 580" role="img" aria-labelledby="title desc">',
719 '<title id="title">Root of unity projection transported into falling factorial coordinates</title>',
720 '<desc id="desc">The fifth Touchard polynomial with coefficients 1,15,25,10,1 is projected by degree modulo three. After falling factorial transport, packet values at x=4 are 600, 244 and 180. They sum exactly to 1024. A degree-five raw-power filter alone gives a different, incorrect packet assignment.</desc>',
721 '<rect width="1080" height="580" fill="#fff"/>',
722 '<text x="40" y="51" font-family="sans-serif" font-size="27" font-weight="700" fill="#16382f">One polynomial, three residue packets, exact recovery</text>',
723 '<text x="40" y="84" font-family="sans-serif" font-size="17" fill="#52665e">T5(z) = z + 15z^2 + 25z^3 + 10z^4 + z^5 | project degrees modulo 3</text>',
724 '<rect x="40" y="113" width="1000" height="67" rx="12" fill="#f1f5f4"/>',
725 '<text x="60" y="143" font-family="sans-serif" font-size="17" font-weight="600" fill="#16382f">Linear coordinate transport</text>',
726 '<text x="60" y="164" font-family="sans-serif" font-size="16" fill="#3b5a50">Phi(z^k) = x(x-1)...(x-k+1) | Q_r = Phi P_r Phi inverse</text>',
727 ]
728 for r in range(3):
729 x=40+340*r
730 out.extend([
731 f'<rect x="{x}" y="214" width="320" height="221" rx="14" fill="{fills[r]}" stroke="{strokes[r]}" stroke-width="1.3"/>',
732 f'<text x="{x+20}" y="251" font-family="sans-serif" font-size="18" font-weight="700" fill="{strokes[r]}">Residue {r} (mod 3)</text>',
733 f'<text x="{x+20}" y="290" font-family="sans-serif" font-size="18" fill="#203a32">{pieces[r]}</text>',
734 f'<text x="{x+20}" y="332" font-family="sans-serif" font-size="15" fill="#486259">At x = 4</text>',
735 f'<text x="{x+20}" y="393" font-family="sans-serif" font-size="49" font-weight="700" fill="{strokes[r]}">{values[r]}</text>',
736 ])
737 out.extend([
738 '<rect x="40" y="462" width="1000" height="80" rx="12" fill="#153f34"/>',
739 '<text x="66" y="496" font-family="sans-serif" font-size="17" fill="#deefeb">Sum of all three transported packets</text>',
740 '<text x="66" y="527" font-family="sans-serif" font-size="30" font-weight="700" fill="#fff">600 + 244 + 180 = 1024 = 4^5</text>',
741 '</svg>'
742 ])
743 return "\n".join(out)+"\n"
744
745
746
747 def svg_image_histogram():
748 """Standalone accessible SVG using independently enumerated exact counts."""
749 hist = image_count_histogram(5,4)
750 fills = ("#1d6855","#3b68a2","#b87321")
751 result = [
752 '<svg xmlns="http://www.w3.org/2000/svg" width="1020" height="550" viewBox="0 0 1020 550" role="img" aria-labelledby="title desc">',
753 '<title id="title">Exact image-size histogram for functions from five inputs to four labels</title>',
754 '<desc id="desc">Enumerating all 1024 functions from a five-element set to a four-element set gives 0, 4, 180, 600, 240 and 0 functions using k=0 through k=5 different labels. Grouping k modulo three gives 600, 244 and 180.</desc>',
755 '<rect width="1020" height="550" fill="#fff"/>',
756 '<text x="42" y="49" font-family="sans-serif" font-size="26" font-weight="700" fill="#173e32">Count actual functions by their image size</text>',
757 '<text x="42" y="82" font-family="sans-serif" font-size="16" fill="#48645d">Every map from 5 inputs to 4 labels, grouped by how many distinct labels occur</text>',
758 '<line x1="74" y1="421" x2="949" y2="421" stroke="#637870" stroke-width="2"/>',
759 ]
760 for k,value in enumerate(hist):
761 x=91+k*145
762 height=round(value*274/600)
763 display_height=max(3,height) if value else 0
764 y=421-display_height
765 result.append(
766 f'<rect x="{x}" y="{y}" width="84" height="{display_height}" rx="4" fill="{fills[k%3]}"/>')
767 result.append(
768 f'<text x="{x+42}" y="{max(115,y-12)}" text-anchor="middle" font-family="sans-serif" font-size="21" font-weight="700" fill="#203e34">{value}</text>')
769 result.append(
770 f'<text x="{x+42}" y="452" text-anchor="middle" font-family="sans-serif" font-size="16" fill="#48645d">k = {k}</text>')
771 result.extend([
772 '<text x="43" y="493" font-family="sans-serif" font-size="16" fill="#42675a">Residue 0: 600</text>',
773 '<text x="375" y="493" font-family="sans-serif" font-size="16" fill="#365d91">Residue 1: 244</text>',
774 '<text x="708" y="493" font-family="sans-serif" font-size="16" fill="#955c20">Residue 2: 180</text>',
775 '<text x="43" y="527" font-family="sans-serif" font-size="17" font-weight="700" fill="#173e32">600 + 244 + 180 = 1024 functions</text>',
776 '</svg>'
777 ])
778 return "\n".join(result)+"\n"
779
780
781
782 def svg_translation_carry():
783 """Vector diagram of the cyclic two-tap transfer; values are recomputed."""
784 _, poly, source = example_row()
785 originals=[int(evaluate(q,4)) for q in source]
786 carries=[int(evaluate(delta(source[(r+1)%3]),4)) for r in range(3)]
787 targets=[int(evaluate(project_newton(translate(poly,1),3,r),4))
788 for r in range(3)]
789 naive=[int(evaluate(translate(q,1),4)) for q in source]
790 assert originals==[600,244,180] and carries==[961,240,900]
791 assert targets==[1561,484,1080] and naive==[1500,1205,420]
792 palette=("#16775f","#3c659d","#aa7021")
793 out=[
794 '<svg xmlns="http://www.w3.org/2000/svg" width="1080" height="600" viewBox="0 0 1080 600" role="img" aria-labelledby="title desc">',
795 '<title id="title">Cross-packet carry when translating the fifth power</title>',
796 '<desc id="desc">For three falling-factorial residue packets at x equals four, old values 600,244,180 carry forward difference contributions 961,240,900 from the next residue. The projected shifted polynomial has packet values 1561,484,1080; direct shifts of the old packets instead give 1500,1205,420. Both totals are 3125, but their packet assignments are different.</desc>',
797 '<rect width="1080" height="600" fill="#fff"/>',
798 '<text x="42" y="51" font-family="sans-serif" font-size="27" font-weight="700" fill="#153e33">Translating a polynomial moves information between packets</text>',
799 '<text x="42" y="82" font-family="sans-serif" font-size="17" fill="#4c665c">Q_r T_1 f = Q_r f + Delta Q_(r+1) f (residues mod 3)</text>',
800 '<rect x="40" y="107" width="1000" height="66" rx="10" fill="#edf5f2"/>',
801 '<text x="59" y="136" font-family="sans-serif" font-size="17" fill="#244d3e">Input: f(x) = x^5, evaluated at x = 4. Output: Q_r(f(x+1)) at x = 4.</text>',
802 '<text x="59" y="159" font-family="sans-serif" font-size="16" fill="#44665b">Only the next input residue contributes to each output carry.</text>',
803 ]
804 for r in range(3):
805 xx=40+340*r
806 out.extend([
807 f'<rect x="{xx}" y="197" width="320" height="253" rx="14" fill="#fafcfc" stroke="{palette[r]}" stroke-width="1.5"/>',
808 f'<text x="{xx+17}" y="234" font-family="sans-serif" font-size="19" font-weight="700" fill="{palette[r]}">Target residue {r}</text>',
809 f'<text x="{xx+17}" y="273" font-family="sans-serif" font-size="17" fill="#344e45">Original packet: {originals[r]}</text>',
810 f'<text x="{xx+17}" y="308" font-family="sans-serif" font-size="17" fill="#344e45">Carry from {(r+1)%3}: +{carries[r]}</text>',
811 f'<text x="{xx+17}" y="374" font-family="sans-serif" font-size="44" font-weight="700" fill="{palette[r]}">{targets[r]}</text>',
812 f'<text x="{xx+17}" y="412" font-family="sans-serif" font-size="15" fill="#52685e">Shift old packet: {naive[r]}</text>'
813 ])
814 out.extend([
815 '<rect x="40" y="473" width="1000" height="88" rx="12" fill="#143d32"/>',
816 '<text x="62" y="506" font-family="sans-serif" font-size="18" fill="#e6f5ee">Transferred destination packets</text>',
817 '<text x="62" y="544" font-family="sans-serif" font-size="30" font-weight="700" fill="#fff">1561 + 484 + 1080 = 3125 = 5^5</text>',
818 '</svg>'
819 ])
820 return "\n".join(out)+"\n"
821
822
823 def save(directory: Path, counts):
824 directory.mkdir(parents=True,exist_ok=True)
825 for name in FILES:
826 if (directory/name).exists():
827 raise FileExistsError("Refusing to overwrite "+str(directory/name))
828 second=stirling_second(MAX_DEGREE)
829 with (directory/"newton-stirling.csv").open("w",newline="",encoding="utf-8") as f:
830 w=csv.writer(f,lineterminator="\n")
831 w.writerow(("n","k","stirling_second","newton_delta_k_at_zero_over_k_factorial"))
832 for n in range(MAX_DEGREE+1):
833 c=to_newton(power(n))
834 for k in range(n+1):
835 w.writerow((n,k,second[n][k],str(c[k])))
836 with (directory/"image-histograms.csv").open("w",newline="",encoding="utf-8") as f:
837 w=csv.writer(f,lineterminator="\n")
838 w.writerow(("domain_size","codomain_size","distinct_image_size",
839 "enumerated_functions","stirling_falling_formula"))
840 for codomain in range(5):
841 for n in range(7):
842 hist=image_count_histogram(n,codomain)
843 for k in range(n+1):
844 w.writerow((n,codomain,k,hist[k],
845 str(second[n][k]*evaluate(falling(k),codomain))))
846 with (directory/"packet-egfs.csv").open("w",newline="",encoding="utf-8") as f:
847 w=csv.writer(f,lineterminator="\n")
848 w.writerow(("codomain_parameter","modulus","residue","degree",
849 "ordinary_egf_coefficient","weighted_exact_count"))
850 for r in range(3):
851 series=packet_egf_coefficients(4,3,r,MAX_DEGREE)
852 for n in range(MAX_DEGREE+1):
853 w.writerow((4,3,r,n,str(series[n]),
854 str(series[n]*factorial(n))))
855 with (directory/"spectral-projectors.csv").open("w",newline="",encoding="utf-8") as f:
856 w=csv.writer(f,lineterminator="\n")
857 w.writerow(("power_degree","modulus","residue",
858 "spectral_power_coefficients","conjugated_power_coefficients"))
859 for n in range(7):
860 for r in range(3):
861 spectral=spectral_project_newton(power(n),3,r)
862 conjugated=project_newton(power(n),3,r)
863 w.writerow((n,3,r,";".join(map(str,spectral)),
864 ";".join(map(str,conjugated))))
865 row,p,packets=example_row()
866 with (directory/"translation-blocks.csv").open("w",newline="",encoding="utf-8") as f:
867 w=csv.writer(f,lineterminator="\n")
868 w.writerow(("shift","source_residue","target_residue",
869 "block_polynomial_power_coefficients","block_value_at_4"))
870 for h in (Q(0),Q(1),Q(2)):
871 for source in range(3):
872 for target in range(3):
873 block=translate_packet_block(p,3,source,target,h)
874 w.writerow((str(h),source,target,
875 ";".join(map(str,block)),
876 str(evaluate(block,4))))
877 with (directory/"shifted-packet-ledger.csv").open("w",newline="",encoding="utf-8") as f:
878 w=csv.writer(f,lineterminator="\n")
879 w.writerow(("residue","Qr_f_at_4","T1_Qr_f_at_4",
880 "Delta_Qr_f_at_4","Delta_Q_next_f_at_4",
881 "Qr_T1_f_at_4"))
882 for r in range(3):
883 original=evaluate(packets[r],4)
884 individual=evaluate(translate(packets[r],1),4)
885 carry=evaluate(delta(packets[(r+1)%3]),4)
886 projected=evaluate(project_newton(translate(p,1),3,r),4)
887 self_difference=evaluate(delta(packets[r]),4)
888 w.writerow((r,str(original),str(individual),str(self_difference),
889 str(carry),str(projected)))
890 with (directory/"finite-characteristic-obstruction.csv").open("w",newline="",encoding="utf-8") as f:
891 w=csv.writer(f,lineterminator="\n")
892 w.writerow(("base_field","falling_degree","N_eigenvalue_mod_3",
893 "polynomial_mod_3","even_projector_mod_3",
894 "odd_projector_mod_3"))
895 for k in (0,3):
896 fpoly=falling(k)
897 w.writerow(("F3",k,k%3,";".join(map(str,mod_polynomial(fpoly,3))),
898 ";".join(map(str,mod_polynomial(project_newton(fpoly,2,0),3))),
899 ";".join(map(str,mod_polynomial(project_newton(fpoly,2,1),3)))))
900 with (directory/"packet-contributions.csv").open("w",newline="",encoding="utf-8") as f:
901 w=csv.writer(f,lineterminator="\n")
902 w.writerow(("residue","k","stirling_count","falling_4_k","contribution","character_trace"))
903 for r in range(3):
904 for k,c in enumerate(row):
905 if k%3==r:
906 vk=evaluate(falling(k),4)
907 w.writerow((r,k,c,str(vk),str(c*vk),
908 str(trace_character_selector(3,k,r))))
909 with (directory/"difference-transfer.csv").open("w",newline="",encoding="utf-8") as f:
910 w=csv.writer(f,lineterminator="\n")
911 w.writerow(("residue","packet_value_x4","delta_packet_x4",
912 "projected_delta_in_predecessor_packet_x4"))
913 for r in range(3):
914 w.writerow((r,str(evaluate(packets[r],4)),
915 str(evaluate(delta(packets[r]),4)),
916 str(evaluate(project_newton(delta(p),3,(r-1)%3),4))))
917 with (directory/"refinement.csv").open("w",newline="",encoding="utf-8") as f:
918 w=csv.writer(f,lineterminator="\n")
919 w.writerow(("parent_modulus","parent_residue","child_modulus",
920 "child_residue","child_value_x4","parent_value_x4"))
921 for r in range(3):
922 for q in range(2):
923 w.writerow((3,r,6,r+3*q,
924 str(evaluate(project_newton(p,6,r+3*q),4)),
925 str(evaluate(packets[r],4))))
926 with (directory/"projector-matrix.csv").open("w",newline="",encoding="utf-8") as f:
927 w=csv.writer(f,lineterminator="\n")
928 w.writerow(("residue","row_monomial_degree","column_input_power",
929 "matrix_coefficient"))
930 for r in range(3):
931 for col in range(6):
932 proj=project_newton(power(col),3,r)
933 for row_number in range(6):
934 w.writerow((r,row_number,col,str(
935 proj[row_number] if row_number<len(proj) else Q(0))))
936 (directory/"projection-diagram.svg").write_text(svg_diagram(),encoding="utf-8")
937 (directory/"function-image-histogram.svg").write_text(svg_image_histogram(),encoding="utf-8")
938 (directory/"translation-carry.svg").write_text(svg_translation_carry(),encoding="utf-8")
939 summary={
940 "arithmetic":"exact integer and Fraction",
941 "scope":"Polynomials over Q; finite degree and m>=2; no analytic or p-adic convergence",
942 "row5":list(row),
943 "newton_polynomial":"z + 15z^2 + 25z^3 + 10z^4 + z^5",
944 "packet_polynomials_power_basis":[list(map(str,part)) for part in packets],
945 "packet_values_at_x4":[str(evaluate(v,4)) for v in packets],
946 "image_count_histogram_n5_x4":list(image_count_histogram(5,4)),
947 "number_operator_on_x_squared":list(map(str,newton_degree_operator(power(2)))),
948 "shift_one_project_after_translate_at_4":[str(evaluate(
949 project_newton(translate(p,1),3,r),4)) for r in range(3)],
950 "shift_one_translate_existing_packet_at_4":[str(evaluate(
951 translate(packets[r],1),4)) for r in range(3)],
952 "shift_one_cross_packet_carries_at_4":[str(evaluate(
953 delta(packets[(r+1)%3]),4)) for r in range(3)],
954 "finite_characteristic_obstruction":{
955 "field":"F3","character_modulus":2,
956 "colliding_degree_eigenvalues_mod_3":[0,0],
957 "falling_degree_3_mod_3":list(mod_polynomial(falling(3),3)),
958 "even_projector_on_falling_0_mod_3":[1],
959 "even_projector_on_falling_3_mod_3":[0],
960 },
961 "packet_egf_n5_x4":[str(packet_egf_coefficients(4,3,r,5)[5]*120)
962 for r in range(3)],
963 "difference_values_at_x4":[str(evaluate(delta(v),4)) for v in packets],
964 "reconstructed_value_at_x4":str(evaluate(p,4)),
965 "naive_residue2_at_x4":str(evaluate(naive_power_project(p,3,2),4)),
966 "correct_residue2_at_x4":str(evaluate(packets[2],4)),
967 "star_x_x_coefficients":list(map(str,star(power(1),power(1)))),
968 "checks_by_group":counts,
969 "total_assertions":sum(counts.values())
970 }
971 (directory/"results.json").write_text(json.dumps(summary,indent=2,ensure_ascii=False)+"\n",encoding="utf-8")
972
973
974 def main():
975 parser=argparse.ArgumentParser(description=__doc__)
976 parser.add_argument("--output",type=Path,help="New destination for exact CSV, SVG and JSON")
977 args=parser.parse_args()
978 counts=check()
979 print("Stirling row five:",ROW_FIVE)
980 print("Root-of-unity packets at x=4: 600 + 244 + 180 = 1024")
981 print("Finite-difference packets at x=4: 900 + 961 + 240 = 2101")
982 print("Independent function-image counts: 0, 4, 180, 600, 240, 0")
983 print("After shift, projected packets at x=4: 1561 + 484 + 1080 = 3125")
984 print("Shifted old packets at x=4: 1500 + 1205 + 420 = 3125")
985 print("F3 obstruction: N has a shared eigenvalue on degrees 0 and 3, but Q0 mod 2 separates them")
986 print("Spectral number operator: N=x*(I-backshift), N(x^2)=2x^2-x")
987 print("Transferred projector != raw power-degree filter (180 != 1024)")
988 print("PASS:",sum(counts.values()),"exact assertions")
989 if args.output:
990 save(args.output,counts)
991 print("Exact artifacts:",args.output)
992
993
994 if __name__=="__main__":
995 main()