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