figures_sequences.py
Special Functions
figures_sequences.py
This is the original source text. Reading it does not run or verify the program.
1 """Sixteen exact sequence lessons, with bounded reproducibility evidence.
2
3 The sequences are integer/rational objects. Floating arithmetic is used only to
4 place marks and evaluate explicitly labeled logarithmic display coordinates.
5 This file proves no novelty claim, asymptotic claim, or universal certificate.
6 Run: python figures_sequences.py
7 """
8 from pathlib import Path
9 from fractions import Fraction as Q
10 from math import comb, factorial, gcd, log10
11 from itertools import permutations
12 from datetime import datetime, timezone
13 import hashlib
14 import json
15 import platform
16 import time
17 import numpy as np
18 import matplotlib
19 import matplotlib.pyplot as plt
20 from guide_style import ROOT, COLORS, INK, MUTED, PALE, save_figure
21
22 N = 18
23 STEMS = [
24 'S01_fibonacci_lucas','S02_tribonacci_tetranacci','S03_jacobsthal',
25 'S04_sqrt3_pell','S05_binomial_catalan','S06_catalan_convolution',
26 'S07_bell_stirling','S08_fubini_factorial','S09_derangements',
27 'S10_euler_updown','S11_bernoulli_pair','S12_partition_product',
28 'S13_apery_recurrence','S14_divisor_convolution','S15_binary_families',
29 'S16_polynomial_ladder']
30
31 def recurrence(seed, taps, degree=N):
32 a = list(seed)
33 for n in range(len(a), degree + 1):
34 a.append(sum(t*a[n-j] for j,t in enumerate(taps, 1)))
35 return a
36
37 def convolution(a,b,degree):
38 return [sum(a[k]*b[n-k] for k in range(n+1)) for n in range(degree+1)]
39
40 def invert(a,degree):
41 assert a[0] != 0
42 b = [Q(1,a[0])]
43 for n in range(1,degree+1):
44 b.append(-sum(a[k]*b[n-k] for k in range(1,n+1))/a[0])
45 return b
46
47 def table(ax, rows, labels, cols, font=12, scale=1.7):
48 ax.axis('off')
49 tab=ax.table(cellText=rows,rowLabels=labels if labels else None,colLabels=cols,
50 loc='center',cellLoc='center',rowLoc='center')
51 tab.auto_set_font_size(False);tab.set_fontsize(font);tab.scale(1,scale)
52 for (r,c),cell in tab.get_celld().items():
53 cell.set_edgecolor('#D8E3E4');cell.set_linewidth(.7)
54 if r==0 or c==-1:
55 cell.set_facecolor(PALE);cell.get_text().set_color(COLORS[0])
56 return tab
57
58 def basic_axes(ax, title, xlabel='Index n', ylabel='Value'):
59 ax.set_title(title,pad=10);ax.set_xlabel(xlabel);ax.set_ylabel(ylabel)
60 ax.grid(axis='y');ax.set_axisbelow(True)
61
62 def logvals(a):
63 return np.array([log10(1+x) for x in a])
64
65 def factor_mu(n):
66 original=n; count=0; p=2
67 while p*p<=n:
68 if n%p==0:
69 n//=p;count+=1
70 if n%p==0:return 0
71 while n%p==0:n//=p
72 p+=1
73 if n>1:count+=1
74 return (-1)**count
75
76 def partitions(n,maxpart=None):
77 if n==0:return [()]
78 if maxpart is None:maxpart=n
79 out=[]
80 for j in range(min(n,maxpart),1-1,-1):
81 out.extend([(j,)+p for p in partitions(n-j,j)])
82 return out
83
84 def main():
85 started=datetime.now(timezone.utc).isoformat();t0=time.monotonic()
86 input_names=['figures_sequences.py','chapters_sequences.tex','guide_style.py',
87 'evidence/OEIS_SOURCE_RECORDS.json']
88 inputs=[{'path':p,'sha256':hashlib.sha256((ROOT/p).read_bytes()).hexdigest()}
89 for p in input_names]
90 C={};prefixes={}
91 def record_prefix(identifier, values, first=0):
92 prefixes[identifier]={'initial_index':first,'first_eight_terms':list(values[:8])}
93
94 # 01: one recurrence, two initial states; common-index norm identity.
95 fib=recurrence([0,1],[1,1]);luc=recurrence([2,1],[1,1])
96 defect=[luc[n]**2-5*fib[n]**2 for n in range(N+1)]
97 assert defect==[4*(-1)**n for n in range(N+1)]
98 record_prefix('A000045',fib);record_prefix('A000032',luc)
99 fig,axs=plt.subplots(1,2,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'width_ratios':[1.4,1]})
100 n=np.arange(13)
101 axs[0].plot(n,fib[:13],'-o',label='Fibonacci: seeds 0, 1',ms=4)
102 axs[0].plot(n,luc[:13],'-s',label='Lucas: seeds 2, 1',ms=4)
103 basic_axes(axs[0],'Same recurrence, different seeds');axs[0].set_xticks([0,3,6,9,12]);axs[0].legend(loc='upper left',fontsize=10.5)
104 axs[1].bar(n,defect[:13],color=[COLORS[i%2] for i in n],width=.72)
105 axs[1].axhline(0,color=MUTED,lw=1)
106 basic_axes(axs[1],'Common-index defect',ylabel=r'$L_n^2-5F_n^2$')
107 axs[1].set(yticks=[-4,0,4],ylim=(-5.6,5.6),xticks=[0,3,6,9,12])
108 save_figure(fig,STEMS[0]);C['S01']={'Fibonacci':fib,'Lucas':luc,'defect':defect,'worked_index':7,'worked_defect':-4}
109
110 # 02: exact OEIS seed conventions; no zero is hidden by log scaling.
111 trib=recurrence([0,0,1],[1,1,1]);tetra=recurrence([0,0,0,1],[1,1,1,1])
112 record_prefix('A000073',trib);record_prefix('A000078',tetra)
113 fig,axs=plt.subplots(1,2,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'width_ratios':[1,1.25]})
114 table(axs[0],[trib[:4],tetra[:4]],['Tribonacci','Tetranacci'],['n=0','1','2','3'],font=11.5,scale=2.3)
115 axs[0].set_title('Fix the seeds before iterating',pad=15)
116 axs[1].plot(range(17),logvals(trib[:17]),'-o',ms=4,label='Three-term memory')
117 axs[1].plot(range(17),logvals(tetra[:17]),'-s',ms=4,label='Four-term memory')
118 basic_axes(axs[1],'Growth on a labeled log scale',ylabel=r'$\log_{10}(1+a_n)$');axs[1].set_xticks([0,4,8,12,16]);axs[1].legend(fontsize=10.5)
119 save_figure(fig,STEMS[1]);C['S02']={'tribonacci':trib,'tetranacci':tetra,'plot_transformation':'log10(1+a_n)'}
120
121 # 03: independently compare a closed form with recurrence iteration.
122 jac=recurrence([0,1],[1,2]);closed=[(2**n-(-1)**n)//3 for n in range(N+1)]
123 assert jac==closed
124 record_prefix('A001045',jac);rem=[2**n-3*jac[n] for n in range(N+1)]
125 fig,axs=plt.subplots(1,2,figsize=(8.1,4.1),layout='constrained')
126 n=np.arange(13);axs[0].bar(n,jac[:13],color=COLORS[0])
127 basic_axes(axs[0],'Jacobsthal growth');axs[0].set_xticks([0,3,6,9,12])
128 axs[1].bar(n,rem[:13],color=[COLORS[i%2] for i in n]);axs[1].axhline(0,color=MUTED,lw=1)
129 basic_axes(axs[1],'Remove the growing mode',ylabel=r'$2^n-3J_n$');axs[1].set(yticks=[-1,0,1],ylim=(-1.4,1.4),xticks=[0,3,6,9,12])
130 save_figure(fig,STEMS[2]);C['S03']={'Jacobsthal':jac,'mode_remainder':rem,'closed_form_matches':True}
131
132 # 04: integer multiplication in Z[sqrt(3)], no floating cancellation.
133 x=[1];y=[0]
134 for n in range(N):x.append(2*x[-1]+3*y[-1]);y.append(x[-2]+2*y[-1])
135 assert all(x[n]**2-3*y[n]**2==1 for n in range(N+1))
136 assert x==recurrence([1,2],[4,-1]) and y==recurrence([0,1],[4,-1])
137 record_prefix('A001075',x);record_prefix('A001353',y)
138 fig,axs=plt.subplots(1,2,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'width_ratios':[1.25,1]})
139 table(axs[0],[[n,x[n],y[n],1] for n in range(8)],[],['n','x(n)','y(n)','norm'],font=11.5,scale=1.5)
140 axs[0].set_title('One unit, two integer coordinates')
141 axs[1].axis('off')
142 axs[1].text(.5,.78,'x(n+1) = 2x(n) + 3y(n)\ny(n+1) = x(n) + 2y(n)',ha='center',va='center',fontsize=14)
143 axs[1].text(.5,.48,'(7, 4) → (26, 15)',ha='center',fontsize=15,color=COLORS[0])
144 axs[1].text(.5,.24,'determinant = 1\nnorm x² − 3y² = 1',ha='center',va='center',fontsize=14)
145 save_figure(fig,STEMS[3]);C['S04']={'x':x,'y':y,'norms':[1]*(N+1),'update_matrix':[[2,3],[1,2]]}
146
147 # 05-06: direct binomial Catalan values and independent convolution.
148 central=[comb(2*n,n) for n in range(N+1)];cat=[central[n]//(n+1) for n in range(N+1)]
149 assert all(central[n]%(n+1)==0 for n in range(N+1))
150 assert all(cat[n+1]==sum(cat[k]*cat[n-k] for k in range(n+1)) for n in range(N))
151 record_prefix('A000984',central);record_prefix('A000108',cat)
152 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'height_ratios':[.8,1.3]})
153 table(axs[0],[central[:7],list(range(1,8)),cat[:7]],['Binomial','Divide by','Catalan'],list(range(7)),font=11.5,scale=1.25)
154 n=np.arange(15);axs[1].plot(n,[log10(v) for v in central[:15]],'-o',ms=4,label='Central binomial')
155 axs[1].plot(n,[log10(v) for v in cat[:15]],'-s',ms=4,label='Catalan')
156 basic_axes(axs[1],'Normalize by the changing factor n + 1',ylabel=r'$\log_{10}(\mathrm{value})$');axs[1].set_xticks([0,3,6,9,12,14]);axs[1].legend(fontsize=10.5)
157 save_figure(fig,STEMS[4]);C['S05']={'central_binomial':central,'Catalan':cat,'divisors':list(range(1,N+2))}
158 products=[cat[k]*cat[4-k] for k in range(5)]
159 assert products==[14,5,4,5,14] and sum(products)==42
160 fig,ax=plt.subplots(figsize=(8.1,4.1),layout='constrained')
161 ax.bar(range(5),products,color=COLORS[:5],width=.62)
162 for k,v in enumerate(products):ax.text(k,v+.6,str(v),ha='center',fontsize=15)
163 ax.set(xticks=range(5),xticklabels=['1 × 14','1 × 5','2 × 2','5 × 1','14 × 1'],ylim=(0,18.5))
164 basic_axes(ax,'Five products add to C(5) = 42',xlabel=r'Product $C_kC_{4-k}$',ylabel='Contribution')
165 save_figure(fig,STEMS[5]);C['S06']={'worked_target_index':5,'products':products,'sum':42,'convolution_checked_through_target':N}
166
167 # 07-08: exact Stirling recurrence and two independent row-sum recurrences.
168 S=[[0]*(N+1) for _ in range(N+1)];S[0][0]=1
169 for n in range(1,N+1):
170 for k in range(1,n+1):S[n][k]=S[n-1][k-1]+k*S[n-1][k]
171 bell=[sum(r) for r in S];fub=[sum(factorial(k)*S[n][k] for k in range(n+1)) for n in range(N+1)]
172 bell_alt=[1];fub_alt=[1]
173 for n in range(1,N+1):
174 bell_alt.append(sum(comb(n-1,k)*bell_alt[k] for k in range(n)))
175 fub_alt.append(sum(comb(n,j)*fub_alt[n-j] for j in range(1,n+1)))
176 assert bell==bell_alt and fub==fub_alt
177 record_prefix('A000110',bell);record_prefix('A000670',fub);record_prefix('A000142',[factorial(n) for n in range(N+1)])
178 # Exact EGF checks: exp(e^t-1), via derivative recurrence; reciprocal 2-e^t.
179 h=[Q(0)]+[Q(1,factorial(n)) for n in range(1,N+1)]
180 exp_h=[Q(1)]
181 for n in range(1,N+1):exp_h.append(sum(k*h[k]*exp_h[n-k] for k in range(1,n+1))/n)
182 assert exp_h==[Q(bell[n],factorial(n)) for n in range(N+1)]
183 denominator=[Q(1)]+[-Q(1,factorial(n)) for n in range(1,N+1)]
184 assert invert(denominator,N)==[Q(fub[n],factorial(n)) for n in range(N+1)]
185 fig,axs=plt.subplots(1,2,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'width_ratios':[3.2,1]})
186 A=np.array([[S[n][k] if k<=n else np.nan for k in range(8)] for n in range(8)],float)
187 axs[0].imshow(np.log1p(A),cmap='Blues',aspect='auto')
188 for n in range(8):
189 for k in range(n+1):axs[0].text(k,n,str(S[n][k]),ha='center',va='center',fontsize=11,color='white' if S[n][k]>25 else INK)
190 axs[0].set(xticks=range(8),yticks=range(8),xlabel='Number of blocks k',ylabel='Set size n',title='Stirling triangle S(n, k)')
191 axs[1].axis('off');axs[1].set(xlim=(0,1),ylim=(7.5,-.5))
192 for n in range(8):
193 axs[1].text(.5,n,str(bell[n]),ha='center',va='center',fontsize=13,color=COLORS[0],bbox={'facecolor':PALE,'edgecolor':'none','boxstyle':'round,pad=.28'})
194 axs[1].set_title('Bell row sum B(n)')
195 save_figure(fig,STEMS[6]);C['S07']={'Stirling_rows':S,'Bell':bell,'independent_recurrence_match':True,'EGF_coefficient_match':True,'shade':'log1p; numbers exact'}
196 weights=[factorial(k) for k in range(1,6)];stir=S[5][1:6];weighted=[a*b for a,b in zip(stir,weights)]
197 assert weighted==[1,30,150,240,120] and sum(weighted)==541
198 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'height_ratios':[.8,1.2]})
199 table(axs[0],[stir,weights,weighted],['S(5, k)','k!','Product'],list(range(1,6)),font=12,scale=1.25)
200 axs[1].bar(range(1,6),weighted,color=COLORS[:5],width=.6)
201 for k,v in enumerate(weighted,1):axs[1].text(k,v+6,str(v),ha='center',fontsize=12)
202 basic_axes(axs[1],'Ordered blocks: contributions sum to 541',xlabel='Number of blocks k',ylabel='Weighted count');axs[1].set(xticks=range(1,6),ylim=(0,285))
203 save_figure(fig,STEMS[7]);C['S08']={'Fubini':fub,'factorial_weights':weights,'worked_products':weighted,'sum':541,'independent_recurrence_match':True,'EGF_coefficient_match':True}
204
205 # 09: signed sum, recurrence, and brute-force counts in a deliberately small range.
206 der=[sum((-1)**k*(factorial(n)//factorial(k)) for k in range(n+1)) for n in range(N+1)]
207 alt=[1,0]
208 for n in range(2,N+1):alt.append((n-1)*(alt[n-1]+alt[n-2]))
209 assert der==alt
210 brute=[sum(all(i!=p[i] for i in range(n)) for p in permutations(range(n))) for n in range(9)]
211 assert brute==der[:9]
212 record_prefix('A000166',der)
213 terms=[(-1)**k*(factorial(6)//factorial(k)) for k in range(7)]
214 fig,ax=plt.subplots(figsize=(8.1,4.1),layout='constrained')
215 ax.bar(range(7),terms,color=[COLORS[0] if v>0 else COLORS[1] for v in terms],width=.62)
216 for k,v in enumerate(terms):ax.text(k,v+(35 if v>0 else -40),str(v),ha='center',va='bottom' if v>0 else 'top',fontsize=12)
217 ax.axhline(0,color=MUTED,lw=1);ax.axhline(265,color=COLORS[2],ls='--',label='Final sum = 265')
218 basic_axes(ax,'Inclusion-exclusion at n = 6',xlabel='Excluded fixed-point count k',ylabel='Signed contribution');ax.set(xticks=range(7),ylim=(-910,950));ax.legend(loc='upper right',fontsize=11)
219 save_figure(fig,STEMS[8]);C['S09']={'derangements':der,'terms_at_6':terms,'brute_force_max_n':8,'brute_force_counts':brute,'recurrence_domain':'n>=2'}
220
221 # 10: formal rational inversion of cosine; sine times secant.
222 cosine=[Q(0)]*(N+1);sine=[Q(0)]*(N+1)
223 for n in range(N+1):
224 if n%2==0:cosine[n]=Q((-1)**(n//2),factorial(n))
225 else:sine[n]=Q((-1)**((n-1)//2),factorial(n))
226 secant=invert(cosine,N);tangent=convolution(sine,secant,N)
227 euler=[int((secant[n]+tangent[n])*factorial(n)) for n in range(N+1)]
228 assert all((secant[n]+tangent[n])*factorial(n)==euler[n] for n in range(N+1))
229 euler_alt=[1,1]
230 for n in range(1,N):
231 value=sum(comb(n,k)*euler_alt[k]*euler_alt[n-k] for k in range(n+1))
232 assert value%2==0;euler_alt.append(value//2)
233 assert euler==euler_alt
234 updown_brute=[sum(all((p[i]<p[i+1]) if i%2==0 else (p[i]>p[i+1]) for i in range(n-1)) for p in permutations(range(n))) for n in range(9)]
235 assert updown_brute==euler[:9]
236 record_prefix('A000111',euler);record_prefix('A000364',euler[0::2]);record_prefix('A000182',euler[1::2],1)
237 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained')
238 table(axs[0],[[euler[2*m] for m in range(6)],[factorial(2*m) for m in range(6)]],['Secant number','Divide by'],[f'deg {2*m}' for m in range(6)],font=11.5,scale=1.55)
239 axs[0].set_title('Even degrees belong to secant')
240 table(axs[1],[[euler[2*m+1] for m in range(6)],[factorial(2*m+1) for m in range(6)]],['Tangent number','Divide by'],[f'deg {2*m+1}' for m in range(6)],font=11.5,scale=1.55)
241 axs[1].set_title('Odd degrees belong to tangent')
242 save_figure(fig,STEMS[9]);C['S10']={'Euler_up_down':euler,'secant_numbers':euler[0::2],'tangent_numbers':euler[1::2],'EGF_series_match':True,'independent_differential_recurrence_match':True,'permutation_max_n':8,'permutation_counts':updown_brute}
243
244 # 11: B1 = -1/2. Rational recurrence plus exact EGF reciprocal product.
245 bern=[Q(1)]
246 for n in range(1,N+1):bern.append(-sum(comb(n+1,k)*bern[k] for k in range(n))/Q(n+1))
247 assert bern[1]==-Q(1,2)
248 series=[bern[n]/factorial(n) for n in range(N+1)]
249 g=[Q(1,factorial(n+1)) for n in range(N+1)]
250 assert convolution(series,g,N)==[Q(1)]+[Q(0)]*N
251 numer=[b.numerator for b in bern];denom=[b.denominator for b in bern]
252 record_prefix('A027641',numer);record_prefix('A027642',denom)
253 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained',sharex=True)
254 n=np.arange(13)
255 axs[0].scatter(n,numer[:13],c=[COLORS[0] if v>=0 else COLORS[1] for v in numer[:13]],s=48,zorder=3)
256 axs[0].vlines(n,0,numer[:13],colors=[COLORS[0] if v>=0 else COLORS[1] for v in numer[:13]])
257 axs[0].set_yscale('symlog',linthresh=1);axs[0].set_yticks([-100,-1,0,1,5]);axs[0].set_yticklabels(['−100','−1','0','1','5'])
258 basic_axes(axs[0],'Reduced numerator p(n): signed log scale',xlabel='',ylabel='Numerator');axs[0].set_ylim(-1600,11)
259 axs[0].annotate('−691',xy=(12,-691),xytext=(10.7,-70),fontsize=11.5,ha='center',arrowprops={'arrowstyle':'->','color':MUTED,'lw':1})
260 axs[1].bar(n,denom[:13],color=COLORS[2],width=.58);axs[1].set_yscale('log')
261 basic_axes(axs[1],'Positive denominator q(n): log scale',ylabel='Denominator');axs[1].set_xticks(range(13));axs[1].set_ylim(.8,6500)
262 axs[1].text(12,2730,'2730',ha='right',va='bottom',fontsize=10.5)
263 save_figure(fig,STEMS[10]);C['S11']={'Bernoulli_rationals':list(map(str,bern)),'numerators':numer,'denominators':denom,'B1':'-1/2','exact_EGF_reciprocal_match':True}
264
265 # 12: finite product truncated by degree, independently enumerated small counts.
266 part=[0]*(N+1);part[0]=1
267 for j in range(1,N+1):
268 for n in range(j,N+1):part[n]+=part[n-j]
269 enum=[len(partitions(n)) for n in range(11)]
270 assert part[:11]==enum
271 record_prefix('A000041',part)
272 parts4=partitions(4);assert len(parts4)==5
273 fig,axs=plt.subplots(1,2,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'width_ratios':[.85,1.5]})
274 table(axs[0],[[' + '.join(map(str,p))] for p in parts4],[],['Partitions of 4'],font=13,scale=2)
275 axs[0].set_title('Five choices of multiplicities')
276 axs[1].bar(range(19),part,color=COLORS[0],width=.68)
277 basic_axes(axs[1],'Prefix of the reciprocal product',ylabel='Partition count p(n)');axs[1].set_xticks([0,3,6,9,12,15,18])
278 save_figure(fig,STEMS[11]);C['S12']={'partitions':part,'partitions_of_4':[list(p) for p in parts4],'independent_enumeration_through':10,'enumerated_counts':enum,'finite_factors':list(range(1,N+1))}
279
280 # 13: binomial sum and polynomial-coefficient recurrence are separate code paths.
281 apery=[sum(comb(n,k)**2*comb(n+k,k)**2 for k in range(n+1)) for n in range(N+1)]
282 for n in range(1,N):
283 assert (n+1)**3*apery[n+1]==(2*n+1)*(17*n*n+17*n+5)*apery[n]-n**3*apery[n-1]
284 record_prefix('A005259',apery)
285 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'height_ratios':[1.4,.7]})
286 axs[0].plot(range(15),[log10(v) for v in apery[:15]],'-o',color=COLORS[0],ms=4)
287 basic_axes(axs[0],'Apéry growth: exact integers, logarithmic display',ylabel=r'$\log_{10}A_n$');axs[0].set_xticks([0,3,6,9,12,14])
288 axs[0].text(.02,.96,'A(0..3) = 1, 5, 73, 1445',transform=axs[0].transAxes,ha='left',va='top',fontsize=12)
289 axs[1].axis('off')
290 axs[1].text(.5,.68,'n = 2: 535 × 73 − 8 × 5 = 39015',ha='center',va='center',fontsize=14,color=COLORS[0])
291 axs[1].text(.5,.25,'39015 ÷ 27 = 1445',ha='center',va='center',fontsize=14)
292 save_figure(fig,STEMS[12]);C['S13']={'Apery':apery,'recurrence_n_range':[1,N-1],'worked_n':2,'worked_numerator':39015,'worked_divisor':27,'worked_A3':1445}
293
294 # 14: multiplicative definitions and independent divisor inversion.
295 mu=[0]+[factor_mu(n) for n in range(1,61)]
296 phi=[0]+[sum(gcd(n,k)==1 for k in range(1,n+1)) for n in range(1,61)]
297 mu_alt=[0,1]
298 for n in range(2,61):mu_alt.append(-sum(mu_alt[d] for d in range(1,n) if n%d==0))
299 assert mu==mu_alt
300 assert all(sum(mu[d] for d in range(1,n+1) if n%d==0)==(n==1) for n in range(1,61))
301 assert all(sum(phi[d] for d in range(1,n+1) if n%d==0)==n for n in range(1,61))
302 record_prefix('A008683',mu[1:],1);record_prefix('A000010',phi[1:],1)
303 divs=[d for d in range(1,13) if 12%d==0];mv=[mu[d] for d in divs];pv=[phi[d] for d in divs]
304 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained',sharex=True)
305 axs[0].bar(range(6),mv,color=[COLORS[0] if v>=0 else COLORS[1] for v in mv],width=.6);axs[0].axhline(0,color=MUTED,lw=1)
306 basic_axes(axs[0],'Möbius divisor sum at 12 is zero',xlabel='',ylabel='mu(d)');axs[0].set(yticks=[-1,0,1],ylim=(-1.5,1.5))
307 axs[1].bar(range(6),pv,color=COLORS[2],width=.6)
308 for k,v in enumerate(pv):axs[1].text(k,v+.1,str(v),ha='center',fontsize=12)
309 basic_axes(axs[1],'Totient divisor sum at 12 is twelve',xlabel='Divisor d of 12',ylabel='phi(d)');axs[1].set(xticks=range(6),xticklabels=divs,ylim=(0,5))
310 save_figure(fig,STEMS[13]);C['S14']={'mu_indices_1_60':mu[1:],'phi_indices_1_60':phi[1:],'worked_divisors':divs,'worked_mu':mv,'worked_phi':pv,'mu_sum':sum(mv),'phi_sum':sum(pv),'prime_definition_vs_divisor_inversion_match':True}
311
312 # 15: formulas, channel recurrences, and odd-binomial count for Gould.
313 w=[n.bit_count() for n in range(64)];thue=[v%2 for v in w];gould=[2**v for v in w]
314 for n in range(32):
315 assert w[2*n]==w[n] and w[2*n+1]==w[n]+1
316 assert thue[2*n]==thue[n] and thue[2*n+1]==1-thue[n]
317 assert gould[2*n]==gould[n] and gould[2*n+1]==2*gould[n]
318 odd_count=[sum(comb(n,k)%2 for k in range(n+1)) for n in range(64)]
319 assert gould==odd_count
320 record_prefix('A000120',w);record_prefix('A010060',thue);record_prefix('A001316',gould)
321 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained')
322 for panel,start in enumerate([0,8]):
323 indices=list(range(start,start+8))
324 rows=[[f'{n:04b}' for n in indices],w[start:start+8],thue[start:start+8],gould[start:start+8]]
325 tab=table(axs[panel],rows,['Binary','Weight','Parity','Gould'],indices,font=12,scale=1.6)
326 axs[panel].set_title(f'Binary addresses {start} through {start+7}',pad=13)
327 # Shade odd parity to make the channel pattern legible without extra density.
328 for column,n in enumerate(indices):
329 if thue[n]:tab[(3,column)].set_facecolor('#F4DDD3')
330 save_figure(fig,STEMS[14]);C['S15']={'binary_weight':w,'Thue_Morse':thue,'Gould':gould,'odd_binomial_count_match':True,'channel_recurrences_parent_range':[0,31]}
331
332 # 16: boundary-aware prefix sums and backward differences.
333 ints=list(range(N+1));tri=[n*(n+1)//2 for n in ints];tet=[comb(n+2,3) for n in ints];sq=[n*n for n in ints];cube=[n**3 for n in ints]
334 assert all(tri[n]==sum(range(1,n+1)) and tet[n]==sum(tri[1:n+1]) for n in ints)
335 assert all(tri[n]-tri[n-1]==n and tet[n]-tet[n-1]==tri[n] for n in range(1,N+1))
336 assert all(sq[n]-sq[n-1]==2*n-1 and cube[n]-cube[n-1]==3*n*n-3*n+1 for n in range(1,N+1))
337 record_prefix('A000027',ints[1:],1);record_prefix('A000217',tri);record_prefix('A000292',tet);record_prefix('A000290',sq);record_prefix('A000578',cube)
338 fig,axs=plt.subplots(2,1,figsize=(8.1,4.1),layout='constrained',gridspec_kw={'height_ratios':[1.5,.7]})
339 table(axs[0],[ints[1:7],tri[1:7],tet[1:7],sq[1:7],cube[1:7]],['n','Triangular','Tetrahedral','Square','Cube'],list(range(1,7)),font=12,scale=1.45)
340 axs[0].set_title('Common indices: n = 1 through 6')
341 axs[1].axis('off')
342 for ypos,color,direction,operation in [(.75,COLORS[0],1,'prefix sum'),(.2,COLORS[1],-1,'difference')]:
343 for xpos,label in [(.08,'n'),(.5,'T(n)'),(.92,'H(n)')]:
344 axs[1].text(xpos,ypos,label,ha='center',va='center',fontsize=14,color=color)
345 for left,right in [(.13,.44),(.56,.86)]:
346 xa,xb=(left,right) if direction==1 else (right,left)
347 axs[1].annotate('',xy=(xb,ypos),xytext=(xa,ypos),arrowprops={'arrowstyle':'->','lw':1.5,'color':color})
348 axs[1].text((left+right)/2,ypos+.12,operation,ha='center',fontsize=10.5,color=color)
349 save_figure(fig,STEMS[15]);C['S16']={'positive_integer_extension_at_zero':0,'triangular':tri,'tetrahedral':tet,'squares':sq,'cubes':cube,'difference_domain':'1<=n<=18','boundary_T0_H0':[0,0]}
350
351 # Exact generated prefix comparison with fresh authenticated OEIS metadata.
352 current=json.loads((ROOT/'evidence/OEIS_SOURCE_RECORDS.json').read_text())
353 meta={r['id']:r for r in current['records']}
354 assert len(prefixes)==30 and set(prefixes)==set(meta)
355 for identifier,p in prefixes.items():
356 assert p['initial_index']==meta[identifier]['initial_index'],identifier+' offset'
357 assert p['first_eight_terms']==meta[identifier]['first_eight_terms'],identifier+' prefix'
358 for item in inputs:item['sha256_after']=hashlib.sha256((ROOT/item['path']).read_bytes()).hexdigest();assert item['sha256']==item['sha256_after']
359 outputs=[]
360 for stem in STEMS:
361 for suffix in ['.pdf','.png']:
362 p=ROOT/'figures'/(stem+suffix)
363 outputs.append({'path':str(p.relative_to(ROOT)),'sha256':hashlib.sha256(p.read_bytes()).hexdigest()})
364 manifest={
365 'schema_version':1,'claim_id':'sequence-atlas-exact-definitions-v0.3',
366 'repository':{'commit':None,'dirty':True,'artifact_version':'v0.3; local document source; no Git revision'},
367 'command':['python','figures_sequences.py'],
368 'environment':{'software':[f'Python {platform.python_version()}',f'NumPy {np.__version__}',f'Matplotlib {matplotlib.__version__}'],'hardware':platform.platform()},
369 'mathematics':{'assertion_tested':'All sixteen displayed finite constructions and worked examples reproduce direct definitions and stated independent checks; all thirty generated OEIS prefixes match pinned official metadata.',
370 'coefficient_domain':'exact Python integers and fractions.Fraction rationals; floating point only for plotting coordinates',
371 'conventions':'OEIS native offsets; B1=-1/2; ordinary versus exponential coefficients stated per lesson; causal/backward differences used only in stated ranges',
372 'inputs':[p['path'] for p in inputs],
373 'bounds':{'default_max_index':18,'sequence_prefix_terms':8,'permutation_max_n':8,'partition_enumeration_max_n':10,'divisor_max_n':60,'binary_max_n':63,'OEIS_source_commit':current['source_commit']},
374 'non_claims':['Finite checks do not substitute for universal proofs of the classical identities.','No claim of novelty, asymptotic growth certification, generic fixed taps for triangular transforms, or Lean formalization.','Plot placement and logarithmic ordinates are approximate display coordinates; sequence arithmetic is exact.']},
375 'randomness':{'used':False,'generator':'none','seed':None},
376 'run':{'started_at':started,'runtime_seconds':time.monotonic()-t0,'exit_status':0},
377 'outputs':outputs,'input_artifacts':inputs,
378 'checks':[{'id':k,'status':'PASS','details':v} for k,v in C.items()]+[{'id':'official_OEIS_prefixes','status':'PASS','count':30}],
379 'sequence_prefixes':prefixes,
380 'result':'implementation and finite assertion verified in the stated range',
381 'residual_risks':['Schema v1 is used for legacy compatibility; explicit input hashes are included as a supplementary field, but no machine-enforced resource cap is claimed.','Final assembled-PDF typography is checked by the parent guide workflow after chapter integration.']}
382 (ROOT/'evidence/sequence_checks.json').write_text(json.dumps(manifest,indent=2)+'\n')
383 print(json.dumps({'sequence_figures':16,'exact_checks':'PASS','official_prefixes':30,'runtime_seconds':manifest['run']['runtime_seconds']}))
384
385 if __name__=='__main__':main()
Other PDFs for this work
Original file and source record
This copy matches the recorded source file. Its publication here does not change the work’s review status.
SHA-256: 9381f0cd14c9ff99a66a2b63b1951e2cfed27461c87c55b7ad4a95f1e1ea0b53