分割数列を串刺しする合同式
このページは[url=https://www.geogebra.org/m/ddm5798j]マス旅[/url]の一部です。[br][br]「整数の分割数」列といえば、オイラーの母関数の発見と説明が有名ですね。[br][br]それで終わりではありませんでした。[br]ラマヌジャンはその分割数列を串刺しにする合同式をみつけたり、[br]分割数列そのものをパラメータを使った関数にしたりと、さらに深堀しました。[br][br]今回は「分割数列」です。
ラマヌジャンに入る前に[br]かるくオイラーの「分割数列」とその説明を確認しておこう。[br][br]「5の和分解」の問題を解くのに、[br]5=5, 4+1, 3+2, 3+1+1, 2+2+1, 2+1+1+1, 1+1+1+1+1 と相互の脈絡なく、数の大小関係だけで7通りを求めるでしょう。[br]オイラーさんなら、[br]1の倍数、2の倍数、3の倍数、4の倍数、5の倍数の指数の[url=https://www.geogebra.org/m/twxxx3yq#material/h45heu5t]母関数[/url]を作り、[br]その積を求めるでしょう。[br]{1+x^1+x^(1+1)+x^(1+1+1)+x^(1+1+1+1)+x^(1+1+1+1+1)+...}×[br]{1+x^2+x^(2+2)+...}×[br]{1+x^3+x^(3+3)+...}×[br]{1+x^4+....}×[br]{1+x^5+...}×......[br]=1+ax^1+bx^2+cx^3+dx^4+ex^5+.......[br][br]ここで、xの絶対値が1より小さいと仮定して、eが7になればよいですよね。[br][br]x進法だと思って係数を6桁分かける。[br]1の位はすべて1になるように反転してかけます。[br]111111×010101×001001×010001×100001×1.....=...753211だから、[br]1の位から順に1,a,b,c,d,e=[b]1,1,2,3,5,7[/b]となる。[br]だから、e=7となるね。つまり、5次の係数が5の和分解7になるのです。[br]だから、[b][color=#0000ff]1,a.b,c,d,eは0,1,2,3,4,5の和分解を表す[/color][/b]ことになるのですね。[br][br]1+y+y^2+y^3+.....=1/(1-y)とかけること、[br]φ(x)=(1-x)(1-x^2)(1-x^3)(1-x^4).........[br]nの和分解をp(n)とかくことにすると、[br][br]無限級数の積の左辺は[br][b]1/(1-x)*1/(1-x^2)*1/(1-x^3)*.....=1/φ(x)となります。[br]右辺はp(0)+p(1)x^1+p(2)x^2+p(3)x^3+p(4)x^5+p(5)x^5+.........です。[br]1/φ(x)は和分解数列の母だったのですね。[/b][br][br]これで終わりではありません。[br]φ(x)を展開しましょう。[br]φ(x)=1-x^1-x^2+x^5+x^7-x^12-x^15+x^22+x^26+.....。[br]この項の偶数番目の指数をならべてみましょう。1,5,12,22です。[br]これはなんだかわかりすか、三角数、四角数ならぱっと出る人が多いでしょうけど、[br]これは5角数です。[br]5角数の公式は知らなくても、少し書き出せば、[b]f(n)=n*(3n-1)/2[/b]が作れるでしょう。[br]n=1,2,3,4を入れると1,5,12,22となるのです。分割数の影に5角数あり。[br][br]話はこれで終わりません。[br]奇数番目の指数は0,2,7,15,26。これはゴミではないのです。[br]f(0)=0,f(-1)=2,f(-2)=7,f(-3)=15,f(-4)=26。[br][br]これがオイラーの見つけた「[b]分割数の構造化と5角数との表裏一体性[/b]」でした。[br][br]では、これ以上に、ラマヌジャンは何を見つけたというのでしょうか?
ラマヌジャンはオイラーの見つけた分割数列p(n)を串刺しにしました。[br][br]数列p(n)をnとペアでならべてみましょう。[br](n,p(n))のタプルです。[br](1,1),(2,2),(3,3),[u](4,5)[/u],[br](5,7),(6,11),[u](7,15)[/u],[br](8,22),[u](9,30)[/u],[br](10,42),(11,56),(12,77),[br](13,101),[u](14,135)[/u],(15,176),[br](16,231),(17,297),.................[br][url=https://www.geogebra.org/m/twxxx3yq#material/zpvazfm9]p進的な距離[/url]のときもそうですが、大切な視点なのは、素数、差、剰余[br]そして、ラマヌジャンはブロードキャスト・串刺しでまとめてみることです。[br]p(n)が5の倍数になるnを抜き出すと、[u]4[/u],7,[u]9,14[/u]です。 [br][b][color=#0000ff][size=150][size=200]n≡4(mod 5)ならp(n)≡0(mod5)[br][/size][/size][/color][/b]というきれいな決まりが見つかりますね。[br]p(n)が7の倍数になるnを抜き出すと、[u]5[/u],10,11,[u]12[/u]です。[br][b][color=#0000ff][size=200]n≡5(mod 7)ならp(n)≡0(mod7)[br][/size][/color][/b]というきれいな決まりが見つかります。[br]p(n)が11の倍数になるnを抜き出すと、[u]6[/u],8,12,15,16,[u]17[/u]です。[br][b][size=200][color=#0000ff]n≡6(mod 11)ならp(n)≡0(mod11)[/color][/size][br][/b]というきれいな決まりが見つかります。[br][br]ただ、1対1の法則というミクロなバカ真面目な視点ではなく、[br]もっと[b]俯瞰した態度[/b]。[br][br][b]必要条件[/b]だけでいい、[br][br]素数が5,7,11と素数の[b]番号が1ずつ増える[/b]と、[br][b]剰余が[/b]4,5,6と[b]1ずつ増える[/b]という[br][br][b]串刺し的、ブロードキャストな発想[/b]が面白いですよね。[br]
[size=150][b]<[/b][b]分割数列を作ろう>[br][/b][/size][br][color=#0000ff][b][size=150]nの分割数をp[n][/size][/b][/color]としよう。[br][br]階段を上るときに1段のぼりと2段のぼりの2通りがあるとしたら、[br]階段の上り方がフィボナッチ数列になることを思い出そう。[br]ある段の上り方は、1つ前の段と2つ前の段の上り方数の和になる。[br]この漸化式の考え方が動的計画法につながる。[br]2段上りまでできるのが、フィボナッチ、[br]3段上りまでできるのが、トリボナッチだった。[br]n段目までは1段前からn段前からすべて上ることができることを利用しよう。[br]ただし、この発想だけだと、いつもn以下の上り方ができるため、和分解にはならない。[br]たとえば、4=1+1+2、1+2+1、2+1+1のように4段上りには順番がつけられる。[br]そこで、[br][b]最初は各段に1段上りだけでいく、[br]次は各段に2段上りをつけたす。[br]3段上りをつけたす。[br]4段上りをつけたす。。。。。[br][/b]このように、「[b][color=#0000ff]上りのサイズ」は最後が大きくなるようにすればよい[br][/color][/b]ですね。[br]そうすると、1,1,2は[b]1+1+2[/b]に決まる。[br][br]では、具体的なコードにつながるように例を作ろう。[br]最初にすべてのp[n]は0にしておこう。[br]0段上りは1通りで、[b]p[0]=1がフィックス。[/b][br]たとえば、p[5]を求めたいとしよう。[br]1段のぼりのドミノ倒しでp[1]から[5]までがすべて0+1=1通りになる。[b]p[1]=1がフィックス[/b]。[br]2段のぼりのドミノ倒しで、+2をする前の分割数を更新するのが2系統できる。[br] p[2]+=p[0]=>1+1=2, p[3]+=p[1]=>1+1=2,p[4]+=p[2]=>1+2=3,p[5]+=p[3]=>1+2=3。[br][b]p[2]=2がフィックス[/b]。[br]3段上りで3系統のドミノ倒しで、更新[br] P[3]+=p[0]=>2+1=3, p[4]+=p[1]=>3+1=4, p[5]+=p[2]=>3+2=5。[b]p[3]=3でフィックス[/b]。[br]4段上りのドミノ倒しで,p[4]+=p[0]=3+1=4, p[5]+=p[1]=5+1=6。[b]p[4]=4でフィックス[/b]。[br]5段上りのドミノ倒しで,[b]p[5]+=p[0]=6+1=7でフィックス[/b]。[br][br]これをコードにすればよいね。[br][br][color=#9900ff][b][size=150][u]課題:分割数列p(n)を動的計画法でもとめ、ラマヌジャンの法5,7,11の串刺しテストをしよう。[br][/u][/size][/b][/color][br]def generate_partitions(max_n):[br] """オイラーの動的計画法アプローチで高速に分割数列p(n)を生成"""[br] dp = [0] * (max_n + 1)[br] dp[0] = 1[br] for [b]num [/b]in range(1, max_n + 1):[br] for i in range([b]num[/b], max_n + 1):[br][b] dp[i] += dp[i - num][br][/b] return dp[br][br]# ==========================================[br]# 遊びのパラメータ:どこまでの項を観察するか[br]max_check = 100[br]# ==========================================[br][br]p_list = generate_partitions(max_check)[br][br]print(" 【5の串】 n ≡ 4 (mod 5) のとき p(n) を 5 で割った余り")[br]for n in range(max_check + 1):[br] if n % 5 == 4:[br] print(f"p({n:2d}) = {p_list[n]:<15} -> 余り: {p_list[n] % 5}")[br][br]print("\n 【7の串】 n ≡ 5 (mod 7) のとき p(n) を 7 で割った余り")[br]for n in range(max_check + 1):[br] if n % 7 == 5:[br] print(f"p({n:2d}) = {p_list[n]:<15} -> 余り: {p_list[n] % 7}")[br][br]print("\n 【11の串】 n ≡ 6 (mod 11) のとき p(n) を 11 で割った余り")[br]for n in range(max_check + 1):[br] if n % 11 == 6:[br] print(f"p({n:2d}) = {p_list[n]:<15} -> 余り: {p_list[n] % 11}")[br][br][OUT][br] 【5の串】 n ≡ 4 (mod 5) のとき p(n) を 5 で割った余り[br]p( 4) = 5 -> 余り: 0[br]p( 9) = 30 -> 余り: 0[br]p(14) = 135 -> 余り: 0[br]p(19) = 490 -> 余り: 0[br]p(24) = 1575 -> 余り: 0[br]p(29) = 4565 -> 余り: 0[br]p(34) = 12310 -> 余り: 0[br]p(39) = 31185 -> 余り: 0[br]p(44) = 75175 -> 余り: 0[br]p(49) = 173525 -> 余り: 0[br]p(54) = 386155 -> 余り: 0[br]p(59) = 831820 -> 余り: 0[br]p(64) = 1741630 -> 余り: 0[br]p(69) = 3554345 -> 余り: 0[br]p(74) = 7089500 -> 余り: 0[br]p(79) = 13848650 -> 余り: 0[br]p(84) = 26543660 -> 余り: 0[br]p(89) = 49995925 -> 余り: 0[br]p(94) = 92669720 -> 余り: 0[br]p(99) = 169229875 -> 余り: 0[br][br] 【7の串】 n ≡ 5 (mod 7) のとき p(n) を 7 で割った余り[br]p( 5) = 7 -> 余り: 0[br]p(12) = 77 -> 余り: 0[br]p(19) = 490 -> 余り: 0[br]p(26) = 2436 -> 余り: 0[br]p(33) = 10143 -> 余り: 0[br]p(40) = 37338 -> 余り: 0[br]p(47) = 124754 -> 余り: 0[br]p(54) = 386155 -> 余り: 0[br]p(61) = 1121505 -> 余り: 0[br]p(68) = 3087735 -> 余り: 0[br]p(75) = 8118264 -> 余り: 0[br]p(82) = 20506255 -> 余り: 0[br]p(89) = 49995925 -> 余り: 0[br]p(96) = 118114304 -> 余り: 0[br][br] 【11の串】 n ≡ 6 (mod 11) のとき p(n) を 11 で割った余り[br]p( 6) = 11 -> 余り: 0[br]p(17) = 297 -> 余り: 0[br]p(28) = 3718 -> 余り: 0[br]p(39) = 31185 -> 余り: 0[br]p(50) = 204226 -> 余り: 0[br]p(61) = 1121505 -> 余り: 0[br]p(72) = 5392783 -> 余り: 0[br]p(83) = 23338469 -> 余り: 0[br]p(94) = 92669720 -> 余り: 0[br]
[b][size=150]<振り返り>[/size][/b][br][br]ラマヌジャンはこのようなPCを使った力業を使うはずもない。[br]それなのにすごい近似公式をハーディらと作ってしまったらしい。[br][br][br][b]p(n)~1/(4n√3) e^(π√(2n/3))[/b]とか、[br][br][b]p(n)=1/(2π√2) d/dn(e^λn/λn) +O(e^Hn*1/2) , λn=√(n-1/24)[/b]などです。[br][br][b]ラマヌジャン・ハーディの公式[/b]ですね。[br][br]p(200),p(243),p(721),p(14031)の計算などを巡って、[br]マクマホン、レーマー、ラドマッハーらも相互に刺激しあって[br]近似式を計算したり、開発しあったようです。[br][br][color=#9900ff][b][u][size=150]課題:動的計画法のベタな方法のp(n)と近似式との差を比べてみよう。[br][/size][/u][/b][/color][br]import math[br]import matplotlib.pyplot as plt[br][br]def [b]generate_partitions[/b](max_n):[br] """オイラーの動的計画法で厳密な分割数p(n)を生成(真値)"""[br] dp = [0] * (max_n + 1)[br] dp[0] = 1[br] for num in range(1, max_n + 1):[br] for i in range(num, max_n + 1):[br] dp[i] += dp[i - num][br] return dp[br][br]def [b]ramanujan_approx[/b](n):[br] """ラマヌジャン・ハーディの簡易近似式"""[br] if n <= 0:[br] return 0[br] # 分母: 4 * n * sqrt(3)[br] denominator = 4 * n * math.sqrt(3)[br] # 指数部: pi * sqrt(2n / 3)[br] exponent = math.pi * math.sqrt((2 * n) / 3.0)[br] # 一発計算![br] return math.exp(exponent) / denominator[br][br]# ==========================================[br]# 遊びのパラメータ:n=10 から n=500 まで比較してみる[br]max_n = 500[br]# ==========================================[br][br]p_true = [b]generate_partitions[/b](max_n)[br][br]print(f"{'n':>5} | {'真値 p(n)':>25} | {'ラマヌジャン近似値':>18} | {'誤差の割合 (%)':>10}")[br]print("-" * 75)[br][br]# 代表的なnの値をピックアップして比較[br]test_points = [10, 50, 100, 200, 300, 500][br]for n in test_points:[br] true_val = p_true[n][br] approx_val = ramanujan_approx(n)[br] # 誤差の割合 = |真値 - 近似値| / 真値 * 100[br] error_percent = (abs(true_val - approx_val) / true_val) * 100[br] [br] print(f"{n:5d} | {true_val:27d} | {approx_val:25.2f} | {error_percent:9.2f}%")[br][br][OUT][br] n | 真値 p(n) | ラマヌジャン近似値 | 誤差の割合 (%)[br]---------------------------------------------------------------------------[br] 10 | 42 | 48.10 | 14.53%[br] 50 | 204226 | 217590.50 | 6.54%[br] 100 | 190569292 | 199280893.35 | 4.57%[br] 200 | 3972999029388 | 4100251432187.85 | 3.20%[br] 300 | 9253082936723602 | 9494094811674990.00 | 2.60%[br] 500 | 2300165032574323995027 | 2346386625611059167232.00 | 2.01%[br][br][b][u][color=#9900ff][size=150]課題:PDとラマヌジャン・ハーディの近似式を視覚化してくらべよう。[br][/size][/color][/u][/b][br]import math[br]import matplotlib.pyplot as plt[br][br]def generate_partitions(max_n):[br] """オイラーの動的計画法で厳密な分割数p(n)を生成(真値)"""[br] dp = [0] * (max_n + 1)[br] dp[0] = 1[br] for num in range(1, max_n + 1):[br] for i in range(num, max_n + 1):[br] dp[i] += dp[i - num][br] return dp[br][br]def ramanujan_approx(n):[br] """ラマヌジャン・ハーディの簡易近似式"""[br] if n <= 0:[br] return 0[br] denominator = 4 * n * math.sqrt(3)[br] exponent = math.pi * math.sqrt((2 * n) / 3.0)[br] return math.exp(exponent) / denominator[br][br]# ==========================================[br]# 遊びのパラメータ:n=10 から n=500 までをプロット[br]max_n = 500[br]# ==========================================[br][br]# データの準備[br]p_true_all = [b]generate_partitions[/b](max_n)[br]n_list = list(range(10, max_n + 1))[br]true_vals = [p_true_all[n] for n in n_list][br]approx_vals = [ramanujan_approx(n) for n in n_list][br]error_percents = [(abs(true_vals[i] - approx_vals[i]) / true_vals[i]) * 100 for i in range(len(n_list))][br][br]# グラフ描画領域の設定[br]fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))[br][br]# ----------------------------------------------------[br]# 左グラフ:真値 vs 近似値(対数スケール)[br]# ----------------------------------------------------[br]ax1.plot(n_list, true_vals, label="True p(n) (Euler DP)", color="#1f77b4", linewidth=2.5)[br]ax1.plot(n_list, approx_vals, label="Ramanujan Approx", color="#ff7f0e", linestyle="--", linewidth=2)[br]ax1.set_yscale("log") # 縦軸を対数スケールに[br]ax1.set_title("p(n) Growth: True vs Ramanujan (Log Scale)", fontsize=12, fontweight="bold")[br]ax1.set_xlabel("n", fontsize=10)[br]ax1.set_ylabel("p(n) (Log Value)", fontsize=10)[br]ax1.grid(True, which="both", linestyle=":", alpha=0.6)[br]ax1.legend(fontsize=10)[br][br]# ----------------------------------------------------[br]# 右グラフ:誤差の割合(%)の収束[br]# ----------------------------------------------------[br]ax2.plot(n_list, error_percents, color="#d62728", linewidth=2)[br]ax2.set_title("Relative Error Rate (%)", fontsize=12, fontweight="bold")[br]ax2.set_xlabel("n", fontsize=10)[br]ax2.set_ylabel("Error (%)", fontsize=10)[br]ax2.set_ylim(0, 15) # 誤差の範囲を0〜15%で見やすく固定[br]ax2.grid(True, linestyle=":", alpha=0.6)[br][br]# レイアウトを整えて表示[br]plt.tight_layout()[br]plt.show()[br]
神がかってますね。[br][br]こんな不規則にみえて、しかも猛烈に増大する分割数列をきれいに、[br]√3、π、指数関数を使って、高い精度で予測できるなんて、[br]絶句です。[br]オイラーは天才でしたが、ラマヌジャンは神の子ですね。