Stirling, Bell, and Fubini / complete exact Python source
Triangular matrices, signed inverse, Bell and ordered Bell numbers, Touchard and exponential-moment families, formal generating functions, a parity transfer system, coupled residue transfers and an exact recoverability boundary. Python 3.10 or later, standard library only.
Read the complete hand derivation Download Python · Source and results ZIP · Source identities
Generated exact files: stirling-rows.csv · sequence-checks.csv · moment-bridge.csv · reconstruction.csv · signed-transfer.csv · obstruction.csv · reciprocal-packets.csv · stirling-transform.svg · results.json
Run python -B stirling_bell_fubini.py for exact checks; add --output stirling-results for downloadable records.
1 #!/usr/bin/env python3
2 """Stirling, Bell and Fubini numbers: exact triangular and reciprocal transfers.
3
4 The calculations use Python's standard library and exact integer/Fraction
5 arithmetic. The published examples are classical identities, not new theorems.
6 """
7 from __future__ import annotations
8
9 import argparse
10 import csv
11 import json
12 from fractions import Fraction
13 from math import comb, factorial
14 from pathlib import Path
15
16 Q = Fraction
17 DEGREE = 10
18 REFERENCE_ROW = (0, 1, 15, 25, 10, 1)
19 ALTERNATIVE_ROW_4 = (0, 5, 2, 7, 1)
20
21
22 def stirling_second(max_n):
23 rows = [(1,)]
24 for n in range(1, max_n + 1):
25 prev = rows[-1]
26 rows.append(tuple((prev[k-1] if 0 <= k-1 < n else 0)
27 + k*(prev[k] if k < n else 0)
28 for k in range(n+1)))
29 return tuple(rows)
30
31
32 def signed_stirling_first(max_n):
33 """Coefficients of x(x-1)...(x-n+1), in ascending powers of x."""
34 rows = [(1,)]
35 for n in range(1, max_n + 1):
36 prev = rows[-1]
37 rows.append(tuple((prev[k-1] if 0 <= k-1 < n else 0)
38 - (n-1)*(prev[k] if k < n else 0)
39 for k in range(n+1)))
40 return tuple(rows)
41
42
43 def lower_apply(rows, weights):
44 weights = tuple(weights)
45 return tuple(sum(row[k]*weights[k] for k in range(n+1))
46 for n, row in enumerate(rows))
47
48
49 def bell_fubini(rows):
50 length = len(rows)
51 bell = lower_apply(rows, (1,)*length)
52 ordered = lower_apply(rows, tuple(factorial(k) for k in range(length)))
53 return bell, ordered
54
55
56 def series_product(a, b):
57 if len(a) != len(b):
58 raise ValueError("Truncated series require equal declared degree")
59 return tuple(sum((a[j]*b[n-j] for j in range(n+1)), Q(0))
60 for n in range(len(a)))
61
62
63
64 def series_compose(outer, inner):
65 """Exact truncated formal composition; inner must vanish at zero."""
66 if len(outer)!=len(inner) or not inner or inner[0]!=0:
67 raise ValueError("Formal composition requires equal degree and zero inner constant")
68 degree=len(outer)-1
69 result=[Q(0)]*(degree+1)
70 power=(Q(1),)+(Q(0),)*degree
71 for k,coefficient in enumerate(outer):
72 for n in range(degree+1):
73 result[n]+=Q(coefficient)*power[n]
74 power=series_product(power,inner)
75 return tuple(result)
76
77
78 def formal_log_one_plus(degree):
79 """Coefficients of log(1+t), exactly to the requested degree."""
80 return (Q(0),)+tuple(Q((-1)**(j+1),j) for j in range(1,degree+1))
81
82
83 def stirling_transform_egf(degree, weights):
84 """EGF of S*weights = A(exp(t)-1), A(u)=sum weights[k]u**k/k!."""
85 weights=tuple(weights)
86 if len(weights)<degree+1:
87 raise ValueError("Need every input coordinate through declared degree")
88 input_egf=tuple(Q(weights[k],factorial(k)) for k in range(degree+1))
89 expm1=(Q(0),)+tuple(Q(1,factorial(k)) for k in range(1,degree+1))
90 return series_compose(input_egf,expm1)
91
92
93 def series_reciprocal(a):
94 if not a or a[0] != 1:
95 raise ValueError("Reciprocal series requires constant term exactly one")
96 out = [Q(1)]
97 for n in range(1, len(a)):
98 out.append(-sum((a[j]*out[n-j] for j in range(1,n+1)), Q(0)))
99 return tuple(out)
100
101
102 def egf_coefficients(values):
103 return tuple(Q(v, factorial(n)) for n, v in enumerate(values))
104
105
106 def exponential_series_of_expm1(degree, z=1):
107 """Formal exp(z*(exp(t)-1)) in exact rational arithmetic."""
108 z=Q(z)
109 unit = (Q(1),) + (Q(0),)*degree
110 u = (Q(0),) + tuple(Q(1, factorial(j)) for j in range(1, degree+1))
111 result = [Q(0)] * (degree+1)
112 power = unit
113 for k in range(degree+1):
114 for n in range(degree+1):
115 result[n] += power[n]*z**k/factorial(k)
116 power = series_product(power, u)
117 return tuple(result)
118
119
120 def ordered_reciprocal_series(degree):
121 """1 / (2-exp(t)) as a formal power series about t=0."""
122 denominator = (Q(1),) + tuple(-Q(1,factorial(n))
123 for n in range(1,degree+1))
124 return series_reciprocal(denominator)
125
126
127
128 def touchard_value(row, z):
129 """Evaluate the exact row polynomial sum S(n,k)*z**k."""
130 z=Q(z)
131 return sum((value*z**k for k,value in enumerate(row)),Q(0))
132
133
134 def exponential_moment_value(row, scale):
135 """Gamma moment of T_n(scale*u); E[u**k]=k!, u~Exp(1)."""
136 scale=Q(scale)
137 return sum((value*factorial(k)*scale**k for k,value in enumerate(row)),Q(0))
138
139
140 def geometric_falling_moment(order, mean):
141 """Exact falling-factorial moment of a geometric variable with given mean."""
142 if order<0:
143 raise ValueError("Moment order must be nonnegative")
144 mean=Q(mean)
145 if mean<0:
146 raise ValueError("Geometric mean must be nonnegative")
147 q=mean/(1+mean)
148 return Q(factorial(order))*q**order/((1+mean)*(1-q)**(order+1))
149
150
151 def geometric_moment_value(row, mean):
152 """Compute the polynomial's geometric raw moment in the falling basis."""
153 return sum((value*geometric_falling_moment(k,mean)
154 for k,value in enumerate(row)),Q(0))
155
156
157
158 def interpolate_integer_samples(samples):
159 """Recover monomial coefficients from P(0),...,P(n) by Newton differences."""
160 if not samples:
161 raise ValueError("At least one polynomial sample is required")
162 degree=len(samples)-1
163 differences=[Q(v) for v in samples]
164 result=[Q(0)]*(degree+1)
165 falling=[Q(1)]
166 for j in range(degree+1):
167 scale=differences[0]/factorial(j)
168 for k,value in enumerate(falling):
169 result[k]+=scale*value
170 differences=[differences[i+1]-differences[i]
171 for i in range(len(differences)-1)]
172 next_falling=[Q(0)]*(len(falling)+1)
173 for k,value in enumerate(falling):
174 next_falling[k]-=j*value
175 next_falling[k+1]+=value
176 falling=next_falling
177 return tuple(result)
178
179
180 def parameterized_reciprocal_series(degree, scale):
181 """Formal 1/(1-scale*(exp(t)-1)); no floating-point operations."""
182 scale=Q(scale)
183 denominator=(Q(1),)+tuple(-scale*Q(1,factorial(n))
184 for n in range(1,degree+1))
185 return series_reciprocal(denominator)
186
187
188 def packet(values, modulus, residue):
189 if modulus < 1 or not 0 <= residue < modulus:
190 raise ValueError("Invalid residue packet")
191 return tuple(v if n%modulus == residue else Q(0)
192 for n,v in enumerate(values))
193
194
195 def bell_recurrence(degree):
196 out = [1]
197 for n in range(degree):
198 out.append(sum(comb(n,k)*out[k] for k in range(n+1)))
199 return tuple(out)
200
201
202 def fubini_recurrence(degree):
203 out=[1]
204 for n in range(1,degree+1):
205 out.append(sum(comb(n,j)*out[n-j] for j in range(1,n+1)))
206 return tuple(out)
207
208
209 def direct_set_partition_histogram(n):
210 """Independent restricted-growth-string enumeration (small n only)."""
211 if n==0:
212 return (1,)
213 totals=[0]*(n+1)
214 def visit(length, highest):
215 if length==n:
216 totals[highest+1]+=1
217 return
218 for label in range(highest+2):
219 visit(length+1, max(highest,label))
220 visit(0,-1)
221 return tuple(totals)
222
223
224 def surjection_count(n,k):
225 """Independent inclusion-exclusion count of maps onto k labelled blocks."""
226 return sum((-1)**j * comb(k,j) * (k-j)**n for j in range(k+1))
227
228
229 def weighted_sum(row):
230 return sum(factorial(k)*c for k,c in enumerate(row))
231
232
233 def test_cases():
234 counters={}
235 def verify(condition, group, explanation):
236 if not condition:
237 raise AssertionError(f"{group}: {explanation}")
238 counters[group]=counters.get(group,0)+1
239
240 S=stirling_second(DEGREE)
241 s=signed_stirling_first(DEGREE)
242 bell,ordered=bell_fubini(S)
243 verify(S[5] == REFERENCE_ROW,"reference","Stirling row at n=5")
244 verify(tuple(factorial(k)*S[5][k] for k in range(6))
245 == (0,1,30,150,240,120),"reference","Factorial-weighted row")
246 verify(bell[5]==52 and ordered[5]==541,"reference","Two row sums")
247 verify(bell[4]==15 and ordered[4]==75,"reference","n=4 Bell/Fubini values")
248 verify(s[5]==(0,24,-50,35,-10,1),"reference","Signed inverse row")
249 verify(sum(s[5][k]*ordered[k] for k in range(6))==120,
250 "reference","Inverted Fubini family yields 5!")
251 verify(sum(s[5][k]*bell[k] for k in range(6))==1,
252 "reference","Inverted Bell family yields a constant input")
253 verify(bell[:9]==(1,1,2,5,15,52,203,877,4140),
254 "reference","First nine Bell numbers")
255 verify(ordered[:9]==(1,1,3,13,75,541,4683,47293,545835),
256 "reference","First nine ordered Bell numbers")
257 verify(tuple(sum(S[4]) for _ in range(1))==(15,),
258 "reference","Original fourth row count")
259 verify(tuple(S[4])==(0,1,7,6,1),"reference","Original n=4 coordinates")
260 verify(sum(ALTERNATIVE_ROW_4)==bell[4]
261 and weighted_sum(ALTERNATIVE_ROW_4)==ordered[4],
262 "reference","Both observed summaries agree at n=4")
263 verify(ALTERNATIVE_ROW_4!=S[4] and
264 all(v>=0 for v in ALTERNATIVE_ROW_4),
265 "reference","Different nonnegative integral coordinate profiles")
266 verify(ordered[5]==375+130+30+5+1,"reference","Ordered Bell binomial recurrence n=5")
267
268 # Independent finite combinatorial enumeration does not call S recurrence.
269 for n in range(0,9):
270 histogram=direct_set_partition_histogram(n)
271 for k in range(n+1):
272 verify(histogram[k]==S[n][k],
273 "partitions",f"Enumerated partitions n={n}, k={k}")
274 # Labelled surjections: onto maps are not the same as unordered partitions.
275 for n in range(0,8):
276 for k in range(n+1):
277 verify(surjection_count(n,k)==factorial(k)*S[n][k],
278 "surjections",f"Inclusion-exclusion n={n}, k={k}")
279
280 for n in range(0,DEGREE+1):
281 verify(S[n][n]==1 and s[n][n]==1,"triangular","Unit diagonal")
282 verify(sum(s[n][k]*bell[k] for k in range(n+1))==1,
283 "triangular","Bell-to-constant exact inverse")
284 verify(sum(s[n][k]*ordered[k] for k in range(n+1))==factorial(n),
285 "triangular","Fubini-to-factorial exact inverse")
286 for k in range(n+1):
287 verify(sum(s[n][j]*S[j][k] for j in range(k,n+1))
288 == int(n==k),"triangular","Exact product s*S=I")
289 verify(sum(S[n][j]*s[j][k] for j in range(k,n+1))
290 == int(n==k),"triangular","Exact product S*s=I")
291
292 be=egf_coefficients(bell)
293 fe=egf_coefficients(ordered)
294 eb=exponential_series_of_expm1(DEGREE)
295 rf=ordered_reciprocal_series(DEGREE)
296 u=(Q(0),)+tuple(Q(1,factorial(j)) for j in range(1,DEGREE+1))
297 for n in range(DEGREE+1):
298 verify(be[n]==eb[n],"generating","exp(exp(t)-1) coefficient")
299 verify(fe[n]==rf[n],"generating","1/(2-exp(t)) coefficient")
300 verify(bell[n]==bell_recurrence(DEGREE)[n],
301 "generating","Independent Bell recurrence")
302 verify(ordered[n]==fubini_recurrence(DEGREE)[n],
303 "generating","Independent ordered Bell recurrence")
304 for k in range(DEGREE+1):
305 power=(Q(1),)+(Q(0),)*DEGREE
306 for _ in range(k):
307 power=series_product(power,u)
308 for n in range(k,DEGREE+1):
309 verify(power[n]/factorial(k)==Q(S[n][k],factorial(n)),
310 "generating","Fixed k partition EGF")
311
312 # Same Stirling row, two different functionals: evaluation and gamma moment.
313 # Check each parameter through a second formal generating-function route.
314 for z in (Q(0),Q(1,2),Q(1),Q(2)):
315 ts=tuple(touchard_value(row,z) for row in S)
316 fs=tuple(exponential_moment_value(row,z) for row in S)
317 t_egf=exponential_series_of_expm1(DEGREE,z)
318 f_egf=parameterized_reciprocal_series(DEGREE,z)
319 for n in range(DEGREE+1):
320 verify(t_egf[n]==Q(ts[n],factorial(n)),
321 "moments","Touchard EGF at parameter "+str(z))
322 verify(f_egf[n]==Q(fs[n],factorial(n)),
323 "moments","Parameterized Fubini EGF at "+str(z))
324 verify(sum(s[n][j]*fs[j] for j in range(n+1))
325 == factorial(n)*z**n,
326 "moments","Signed inverse recovers scaled gamma moments")
327 verify(geometric_moment_value(S[n],z)==fs[n],
328 "moments","Independent geometric moment of Touchard polynomial")
329 verify(geometric_falling_moment(n,z)==factorial(n)*z**n,
330 "moments","Geometric falling factorial identity")
331 for n in range(DEGREE):
332 verify(ts[n+1] == z*sum(comb(n,j)*ts[j] for j in range(n+1)),
333 "moments","Touchard binomial recurrence")
334 for n in range(1,DEGREE+1):
335 verify(fs[n] == z*sum(comb(n,j)*fs[n-j] for j in range(1,n+1)),
336 "moments","Parameterized reciprocal recurrence")
337 verify(touchard_value(S[5],Q(1,2))==Q(257,32),
338 "reference","Touchard half evaluation")
339 verify(exponential_moment_value(S[5],Q(1,2))==Q(91,2),
340 "reference","Half-scaled gamma moment")
341 verify(touchard_value(S[5],Q(2))==454,
342 "reference","Touchard evaluation at two")
343 verify(exponential_moment_value(S[5],Q(2))==9002,
344 "reference","Double-scaled gamma moment")
345 verify(touchard_value(S[5],1)==bell[5],
346 "reference","Touchard evaluation recovers Bell")
347 verify(exponential_moment_value(S[5],1)==ordered[5],
348 "reference","Gamma moment recovers Fubini")
349 verify(parameterized_reciprocal_series(DEGREE,1)==rf,
350 "reference","Parameter one recovers original reciprocal")
351 verify(touchard_value(S[0],0)==exponential_moment_value(S[0],0)==1,
352 "reference","Empty-set boundary for zero parameter")
353
354 # General coefficient grammar: composition and its formal logarithmic inverse.
355 weight_families=(
356 (1,)*(DEGREE+1),
357 tuple(factorial(k) for k in range(DEGREE+1)),
358 tuple(Q((-1)**k*(k+1),3) for k in range(DEGREE+1)),
359 tuple(Q(k*k+1,k+1) for k in range(DEGREE+1)),
360 )
361 log= formal_log_one_plus(DEGREE)
362 for weights in weight_families:
363 outputs=lower_apply(S,weights)
364 generated=stirling_transform_egf(DEGREE,weights)
365 recovered=series_compose(generated,log)
366 for n in range(DEGREE+1):
367 verify(generated[n]==Q(outputs[n],factorial(n)),
368 "composition","General Stirling formal substitution")
369 verify(recovered[n]==Q(weights[n],factorial(n)),
370 "composition","Formal logarithm recovers original weight")
371
372 # Touchard polynomial calculus and the exact binomial-type addition law.
373 for z in (Q(0),Q(1,2),Q(1),Q(2),Q(-1)):
374 for n in range(DEGREE):
375 derivative=sum(k*S[n][k]*z**(k-1) for k in range(1,n+1))
376 verify(touchard_value(S[n+1],z)
377 ==z*(touchard_value(S[n],z)+derivative),
378 "touchard","Touchard differential recurrence")
379 for x,y in ((Q(0),Q(1)),(Q(1),Q(1)),(Q(1,2),Q(2)),(Q(-1),Q(2))):
380 for n in range(DEGREE+1):
381 verify(touchard_value(S[n],x+y)==sum(
382 comb(n,j)*touchard_value(S[j],x)
383 *touchard_value(S[n-j],y) for j in range(n+1)),
384 "touchard","Touchard binomial addition law")
385
386 # Full parameter polynomial is reversible; one or two scalar summaries are not.
387 for n,row in enumerate(S):
388 touchard_samples=tuple(touchard_value(row,j) for j in range(n+1))
389 ordered_samples=tuple(exponential_moment_value(row,j) for j in range(n+1))
390 verify(interpolate_integer_samples(touchard_samples)==tuple(row),
391 "reconstruction","Touchard interpolation recovers the row")
392 ordered_coefficients=interpolate_integer_samples(ordered_samples)
393 verify(tuple(ordered_coefficients[k]/factorial(k)
394 for k in range(n+1))==tuple(row),
395 "reconstruction","Ordered polynomial recovers factorial-scaled row")
396 original,candidate=S[4],ALTERNATIVE_ROW_4
397 verify(touchard_value(original,2)==94
398 and touchard_value(candidate,2)==90,
399 "reconstruction","Another evaluation separates the two profiles")
400 verify(exponential_moment_value(original,2)==730
401 and exponential_moment_value(candidate,2)==746,
402 "reconstruction","Another moment separates the two profiles")
403 verify(interpolate_integer_samples(tuple(
404 touchard_value(candidate,j) for j in range(5)))==candidate,
405 "reconstruction","Interpolate the alternative candidate itself")
406 verify(interpolate_integer_samples(tuple(
407 exponential_moment_value(candidate,j) for j in range(5)))
408 ==tuple(candidate[k]*factorial(k) for k in range(5)),
409 "reconstruction","Recover candidate factorial-weighted coefficients")
410
411 # The mod-m filters satisfy an entire coupled transfer, not just reconstruction.
412 for modulus in range(2,6):
413 for residue in range(modulus):
414 expected=packet(fe,modulus,residue)
415 generated=[Q(int(residue==0))]+[Q(0)]*DEGREE
416 for tap in range(modulus):
417 contribution=series_product(packet(u,modulus,tap),
418 packet(fe,modulus,(residue-tap)%modulus))
419 generated=[generated[n]+contribution[n]
420 for n in range(DEGREE+1)]
421 for n in range(DEGREE+1):
422 verify(expected[n]==generated[n],"packet_transfer",
423 "Coupled residue transfer at modulus "+str(modulus))
424 for residue,expected in enumerate((30,380,131)):
425 verify(sum(comb(5,j)*ordered[5-j] for j in range(1,6)
426 if j%3==residue)==expected,
427 "packet_transfer","Three-residue tap contributions at n=5")
428 verify(30+380+131==ordered[5],"packet_transfer",
429 "All three taps reconstruct Fubini n=5")
430
431 # Even/odd packets of F and U obey an exact 2x2 convolutional transfer.
432 ue,uo=packet(u,2,0),packet(u,2,1)
433 even,odd=packet(fe,2,0),packet(fe,2,1)
434 lhs_even=series_product(ue,even)
435 cross_even=series_product(uo,odd)
436 lhs_odd=series_product(ue,odd)
437 cross_odd=series_product(uo,even)
438 for n in range(DEGREE+1):
439 verify(even[n] == (int(n==0)+lhs_even[n]+cross_even[n]),
440 "packets","Even reciprocal packet transfer")
441 verify(odd[n] == lhs_odd[n]+cross_odd[n],
442 "packets","Odd reciprocal packet transfer")
443 verify(even[n]+odd[n]==fe[n],"packets","Two packets reconstruct")
444 verify(sum(packet(fe,3,r)[n] for r in range(3))==fe[n],
445 "packets","Three packets reconstruct")
446 verify(sum(comb(5,j)*ordered[5-j] for j in (2,4))==135,
447 "reference","Even tap-to-odd output contributions")
448 verify(sum(comb(5,j)*ordered[5-j] for j in (1,3,5))==406,
449 "reference","Odd tap-to-even output contributions")
450 verify(135+406==ordered[5],"reference","Two blocks yield F5")
451 return counters
452
453
454 def panel_svg():
455 """Two exact triangular charts with cell contents from integer recurrence."""
456 S=stirling_second(7)
457 top=112
458 cell=52
459 lefts=(68,570)
460 pieces=[
461 '<svg xmlns="http://www.w3.org/2000/svg" width="1050" height="650" viewBox="0 0 1050 650" role="img" aria-labelledby="title desc">',
462 '<title id="title">Stirling partitions and factorial-weighted ordered partitions</title>',
463 '<desc id="desc">Two triangular charts indexed by set size n in rows zero through seven and block count k in columns zero through seven. Left cells are Stirling counts S(n,k); right cells are k factorial times S(n,k). Cells outside k less than or equal to n are left blank. Exact integer values are printed in every valid cell.</desc>',
464 '<rect width="1050" height="650" fill="#ffffff"/>',
465 '<text x="42" y="37" font-family="sans-serif" font-size="23" font-weight="600" fill="#173b39">One triangular kernel · two counting families</text>',
466 '<text x="42" y="63" font-family="sans-serif" font-size="14" fill="#475c59">Unordered partitions become ordered partitions by multiplying each k-block coordinate by k!</text>',
467 ]
468 for chart,left in enumerate(lefts):
469 label=("S(n,k): partitions","k! S(n,k): ordered blocks")[chart]
470 pieces.append(f'<text x="{left}" y="95" font-family="sans-serif" font-size="17" font-weight="600" fill="#235a54">{label}</text>')
471 for k in range(8):
472 xx=left+k*cell+cell//2
473 pieces.append(f'<text x="{xx}" y="126" text-anchor="middle" font-family="monospace" font-size="12" fill="#355751">{k}</text>')
474 for n,row in enumerate(S):
475 yy=top+21+n*cell
476 pieces.append(f'<text x="{left-16}" y="{yy+31}" text-anchor="end" font-family="monospace" font-size="12" fill="#355751">{n}</text>')
477 for k,value in enumerate(row):
478 number=value*factorial(k) if chart else value
479 xx=left+k*cell
480 # Integer bit-length yields a deterministic graduated scale without floating-point logs.
481 shade=min(6,(number.bit_length()-1)//3) if number else 0
482 palette=("#edf4f2","#dcece7","#c2e0d6","#9acbb9",
483 "#6caf9c","#438a7c","#216859")
484 fill=palette[shade]
485 pieces.append(f'<rect x="{xx}" y="{yy}" width="{cell-2}" height="{cell-2}" rx="3" fill="{fill}"/>')
486 pieces.append(f'<text x="{xx+cell//2}" y="{yy+29}" text-anchor="middle" font-family="monospace" font-size="12" font-weight="600" fill="#112d29">{number}</text>')
487 pieces.extend([
488 '<text x="69" y="593" font-family="sans-serif" font-size="15" fill="#334e48">Row n=5: 0 + 1 + 15 + 25 + 10 + 1 = 52</text>',
489 '<text x="570" y="593" font-family="sans-serif" font-size="15" fill="#334e48">Row n=5: 0 + 1 + 30 + 150 + 240 + 120 = 541</text>',
490 '<text x="68" y="624" font-family="sans-serif" font-size="12" fill="#60716e">All displayed values are exact integers. White space above the triangles is outside k ≤ n.</text>',
491 '</svg>'
492 ])
493 return "\n".join(pieces)+"\n"
494
495
496 def save(folder, checks):
497 folder.mkdir(parents=True,exist_ok=True)
498 S=stirling_second(DEGREE)
499 first=signed_stirling_first(DEGREE)
500 bell,ordered=bell_fubini(S)
501 with (folder/"stirling-rows.csv").open("w",encoding="utf-8",newline="") as stream:
502 writer=csv.writer(stream,lineterminator="\n")
503 writer.writerow(("n","k","partitions","block_order_factorial","ordered_partitions"))
504 for n,row in enumerate(S):
505 for k,value in enumerate(row):
506 writer.writerow((n,k,value,factorial(k),value*factorial(k)))
507 with (folder/"sequence-checks.csv").open("w",encoding="utf-8",newline="") as stream:
508 writer=csv.writer(stream,lineterminator="\n")
509 writer.writerow(("n","bell","fubini","independent_bell","independent_fubini","stirling_inverse_bell","stirling_inverse_fubini","fubini_phase_mod2","fubini_phase_mod3"))
510 bb=bell_recurrence(DEGREE)
511 ff=fubini_recurrence(DEGREE)
512 for n in range(DEGREE+1):
513 writer.writerow((n,bell[n],ordered[n],bb[n],ff[n],
514 sum(first[n][j]*bell[j] for j in range(n+1)),
515 sum(first[n][j]*ordered[j] for j in range(n+1)),
516 n%2,n%3))
517 with (folder/"moment-bridge.csv").open("w",encoding="utf-8",newline="") as stream:
518 writer=csv.writer(stream,lineterminator="\n")
519 writer.writerow(("parameter","n","touchard_evaluation","exponential_moment",
520 "geometric_moment","touchard_egf_coefficient","moment_egf_coefficient",
521 "signed_inverse_recovered_weight"))
522 for scale in (Q(0),Q(1,2),Q(1),Q(2)):
523 values=tuple(exponential_moment_value(row,scale) for row in S)
524 for n,row in enumerate(S):
525 writer.writerow((str(scale),n,str(touchard_value(row,scale)),
526 str(values[n]),str(geometric_moment_value(row,scale)),
527 str(touchard_value(row,scale)/factorial(n)),
528 str(values[n]/factorial(n)),
529 str(sum(first[n][j]*values[j]
530 for j in range(n+1)))))
531 with (folder/"signed-transfer.csv").open("w",encoding="utf-8",newline="") as stream:
532 writer=csv.writer(stream,lineterminator="\n")
533 writer.writerow(("n","k","stirling_second","signed_stirling_first","inverse_product"))
534 for n in range(DEGREE+1):
535 for k in range(n+1):
536 product_value=sum(first[n][j]*S[j][k] for j in range(k,n+1))
537 writer.writerow((n,k,S[n][k],first[n][k],product_value))
538 witness=S[4]
539 with (folder/"obstruction.csv").open("w",encoding="utf-8",newline="") as stream:
540 writer=csv.writer(stream,lineterminator="\n")
541 writer.writerow(("k","true_stirling_row","alternative_candidate","difference","factorial","weighted_difference"))
542 for k in range(5):
543 delta=ALTERNATIVE_ROW_4[k]-witness[k]
544 writer.writerow((k,witness[k],ALTERNATIVE_ROW_4[k],delta,factorial(k),factorial(k)*delta))
545 original=S[4]
546 with (folder/"reconstruction.csv").open("w",encoding="utf-8",newline="") as stream:
547 writer=csv.writer(stream,lineterminator="\n")
548 writer.writerow(("parameter","true_touchard","candidate_touchard",
549 "true_ordered_polynomial","candidate_ordered_polynomial"))
550 for z in range(5):
551 writer.writerow((z,touchard_value(original,z),
552 touchard_value(ALTERNATIVE_ROW_4,z),
553 exponential_moment_value(original,z),
554 exponential_moment_value(ALTERNATIVE_ROW_4,z)))
555 with (folder/"reciprocal-packets.csv").open("w",encoding="utf-8",newline="") as stream:
556 writer=csv.writer(stream,lineterminator="\n")
557 writer.writerow(("n","bell","fubini","egf_bell","egf_fubini","even_packet","odd_packet","mod3_residue","mod3_packet"))
558 be=egf_coefficients(bell)
559 fe=egf_coefficients(ordered)
560 for n in range(DEGREE+1):
561 writer.writerow((n,bell[n],ordered[n],str(be[n]),str(fe[n]),
562 str(packet(fe,2,0)[n]),str(packet(fe,2,1)[n]),
563 n%3,str(packet(fe,3,n%3)[n])))
564 (folder/"stirling-transform.svg").write_text(panel_svg(),encoding="utf-8",newline="\n")
565 record=dict(
566 case="Stirling rows, Bell/Fubini triangular transfers and reciprocal parity packets",
567 degree=DEGREE,
568 bell5=bell[5],
569 fubini5=ordered[5],
570 row5=list(S[5]),
571 weighted_row5=[factorial(k)*S[5][k] for k in range(6)],
572 moment_bridge_n5={
573 str(scale):dict(touchard=str(touchard_value(S[5],scale)),
574 exponential_moment=str(exponential_moment_value(S[5],scale)))
575 for scale in (Q(0),Q(1,2),Q(1),Q(2))
576 },
577 witness_n4=dict(true_row=list(S[4]),alternative_row=list(ALTERNATIVE_ROW_4),
578 bell=bell[4],fubini=ordered[4],
579 difference=[ALTERNATIVE_ROW_4[k]-S[4][k] for k in range(5)]),
580 packet_split_n5=dict(even_degree_tap_contribution=135,
581 odd_degree_tap_contribution=406),
582 packet_split_n5_mod3=dict(residue_0=30,residue_1=380,residue_2=131),
583 reconstruction_n4=dict(
584 true_ordered_samples=[str(exponential_moment_value(S[4],z))
585 for z in range(5)],
586 candidate_ordered_samples=[str(exponential_moment_value(
587 ALTERNATIVE_ROW_4,z)) for z in range(5)],
588 true_touchard_at_2=94,candidate_touchard_at_2=90,
589 recovered_stirling_row=[int(v) for v in interpolate_integer_samples(tuple(
590 touchard_value(S[4],z) for z in range(5)))]),
591 checks_by_group=checks,
592 total_assertions=sum(checks.values()),
593 scope="Exact classical finite combinatorics and formal-series identities; no new Lean replay or universal recovery theorem."
594 )
595 (folder/"results.json").write_text(json.dumps(record,indent=2,sort_keys=True)+"\n",
596 encoding="utf-8",newline="\n")
597
598
599 def main():
600 parser=argparse.ArgumentParser(description=__doc__)
601 parser.add_argument("--output",type=Path,help="Generate exact CSV, JSON and SVG tables")
602 args=parser.parse_args()
603 checks=test_cases()
604 print("STIRLING / BELL / FUBINI | exact integer and rational arithmetic")
605 print("n=5: S=(0,1,15,25,10,1); Bell=52; Fubini=541")
606 print("Inverse Stirling transform recovers 1 and k! input families")
607 print("General EGF substitution and exact polynomial interpolation checked")
608 print("Three residue tap paths at n=5: 30 + 380 + 131 = 541")
609 print("Touchard T5(2)=454; gamma moment F5(2)=9002")
610 print("Reciprocal parity packets at degree five: 135 + 406 = 541")
611 print("n=4: distinct nonnegative rows can share Bell=15 and Fubini=75")
612 print(f"PASS: {sum(checks.values())} exact assertions")
613 if args.output:
614 save(args.output,checks)
615 print("Exact tables and diagram written to "+str(args.output))
616
617
618 if __name__=="__main__":
619 main()