ベルヌーイ数は、スイスの数学者 ヤコブ・ベルヌーイ(1654〜1705)が見いだしたと言われる。
その著書「推測術」で述べているらしいが、これが刊行されたのが 1713年だから、歿後の刊行ということになる。
今回は、ベルヌーイ数と累乗の和の計算との関係、その計算方法、また、よく知られた級数との関連について触れたい。
プログラミングでの計算も取り上げる。 サンプルのrubyプログラムは bern.zip に同梱してある。
ベルヌーイ数は、連続する整数の累乗の和を計算する中で見出されたと聞く。
s1(n) = 1^1 + 2^1 + 3^1 + …… + n^1
s2(n) = 1^2 + 2^2 + 3^2 + …… + n^2
s3(n) = 1^3 + 2^3 + 3^3 + …… + n^3
上の計算からベルヌーイ数がどのように立ち現れてくるのか理解できなかったが、自分なりの納得の仕方をみつけたので以降で記したい。
まずは自然数の和から始めよう。s1(n) = 1 + 2 + 3 + …… + n
この s1(n) を計算する際の分かりやすいロジックは、n+1 に着目するというものだとおもう。
1〜10 までの和 s1(10) の例でいうと、
1〜10 を逆に並べた 10, 9, 8, …… 1 を真下に置くことを考える。
1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10
10 + 9 + 8 + 7 + 6 + 5 + 4 + 3 + 2 + 1
上の行の一つの項と、その真下の項を足し算すると、どの項も 11 になる。
これが10個並んでいるから 11*10 = 110 だが、これは s1(10) を２倍した値なので、２で割ればよい。
一般化して書くと s1(n) = 1/2n(n+1) = 1/2(n^2 + n) となる。
ちょっと興味深いのは n^2 が出てくることだ。
x,y座標のグラフで考えると、自然数の和 s1(x) は、
(1, 1) (2, 2) (3, 3) …… (n, n) の点を結ぶ棒グラフを描き、
その棒が織りなす階段状の面積を足し合わせたものである。
階段の上辺は、ギザギザしながらも全体としては45度の角度を保ちながら上に上っていく。
一方、y = x も45度の角度で上昇する。階段のぎざぎざではなく斜線。
その下側に位置する面積(三角形)は y = x^2 / 2 で求めることができる。
y = x を積分したものが y = x^2 / 2 である。
積分で算出される三角形の面積よりも階段の面積の方が大きいが、似通った値になりそうな感じがしないだろうか。
話が脱線したが、1〜n の和が n^2 に遺伝子として組み込まれているのではないか？という感触は、頭の片すみに置いておいていただきたい。
今のところ根拠はないが、s2(n) が n^3 に、s3(n) が n^4 に関連ありそうだという憶測である。
☆ 注目の式： s1(n-1) = 1/2(n^2 - n)
自然数の和 s1(n) を２項定理を用いて考える。
２項定理は (a+b)^n を展開するとどうなるかを示すものだが、
ここでは (x+1)^2 を取り上げる。
(x+1)^2 = x^2 + 2x + 1
これを２項定理に従って正確に書くと
(x+1)^2 = c(2,0) * x^2 + c(2,1) * x * 1 + c(2,2) * 1^2
ここで c(n,k) というのは n個から k個を取り出す組み合わせのパターン数のことである。数学の用語では「２項係数」というようだ。
c(2,0) = 1, c(2,1) = 2, c(2,2) = 1 という値になる。
c(2, 0) は、２個から何も取り出さないパターン数のことで、イメージしにくいが １とみなす。
「組み合わせパターン数」はベルヌーイ数を考えるときに重要な概念として現れるが、とりあえず２乗を考える今の段階では、単なる定数としておく。
２項定理を用いて 5^2 を考えると 5^2 = (4+1)^2 = 4^2 + 2*4 + 1 となる。
これは 4^2 に 2*4 + 1 = 9 を足し算すると 5^2 になることを示している。
16 + 9 = 25 なので間違っていない。
では、4^2 はどうかというと、4^2 = (3+1)^2 = 3^2 + 2*3 + 1 である。
まどろっこしい書き方はやめて、1^2 〜 5^2 を形式的に書き並べてみよう。
1^2 = 1
2^2 = 1^2 + 2*1 + 1
3^2 = 2^2 + 2*2 + 1
4^2 = 3^2 + 2*3 + 1
5^2 = 4^2 + 2*4 + 1
上の式で注目してもらいたいのは、右辺の２乗の項ではなく、その右隣だ。
つまり 2*1, 2*2, 2*3, 2*4 である。
これらは自然数 1〜4 を２倍したものとなっている。
この「２倍したもの」の和が 5^2 に含まれていることは分かるだろうか。
私がこれを理解するのに少し時間がかかったが、要するにこういうことである。
2^2 = 1+(2+1) だが、これに 2*2+1 を加算すると 3^2 になる。
更に、それに 2*3+1 を加算すれば 4^2
更に更に、それに 2*4+1 を加算すれば 5^2 である。
なので、下の式が成り立つ。
5^2 = 1 + (2*1+1) + (2*2+1) + (2*3+1) + (2*4+1)
5^2 = 2 * s1(4) + 1 * 5
s1(4) = (5^2 - 5) / 2
これを一般化して書くと s1(n-1) = 1/2(n^2 - n) = 1/2n(n-1)
上の式で n に１を加算してやれば s1(n) = 1/2n(n+1) となり、当然だが、「オーソドックスな方法」と同じ結果になる。
結果は同じであるものの、「２項定理の利用」の考え方は、２乗の和、３乗の和、４乗の和に応用できそうだと感じてもらえるだろうか。
☆ 注目の式： s2(n-1) = 1/3(n^3 - 3*s1(n-1) - n)
s2(n-1) = 1/3(n^3 - 3/2n^2 + 1/2n)
２乗の和を考えるには、２項定理で３乗の数を算出してみればいい。
(x+1)^3 = x^3 + 3x^2 + 3x + 1
上の式を利用して 1^3 〜 5^3 を書き並べると次のとおり。
1^3 = 1
2^3 = 1^3 + 3*1^2 + 3*1 + 1
3^3 = 2^3 + 3*2^2 + 3*2 + 1
4^3 = 3^3 + 3*3^2 + 3*3 + 1
5^3 = 4^3 + 3*4^2 + 3*4 + 1
上の列記から 5^3 には s2(4) と s1(4) の両方が含まれていることが分かる。
s1(4) は、前項でみたとおり 1/2 * (5^2 - 5) で求めることができる。
5^3 = 3*s2(4) + 3*s1(4) + 5*1
3*s2(4) = 5^3 - 3*s1(4) - 5 = 5^3 - 3/2 * (5^2 - 5) - 5
s2(4) = 1/3 * (5^3 - 3/2 * 5^2 + 1/2 * 5)
s2(4) = 1/3 * (125 - 75/2 + 5/2) = 30
ここで s2(n-1) を一般化して書くと下のとおり。
s2(n-1) = 1/3(n^3 - 3/2n^2 + 1/2n)
一般化して「nまでのk乗の和」を考えるのはまだ早計だが、その手がかりを覗いてみよう。
「組み合わせパターン数」の c(3,x) を使って先の式を書き換えると、
s2(n-1) = 1/3*(c(3,0)*n^3 + c(3,1)*(-1/2)*n^2 + c(3,2)*(1/6)*n)
ベルヌーイ数がちょっと顔を覗かせていることに気が付くだろうか。
ベルヌーイ数に番号をつけて B0, B1, B2 と表現すれば、
B0 = 1, B1 = -1/2, B2 = 1/6 である。
これらベルヌーイ数を使って式を書き換えると、
s2(n-1) = 1/3*(c(3,0)*B0*n^3 + c(3,1)*B1*n^2 + c(3,2)*B2*n)
上は、s2(n-1) を a*n^3 + b*n^2 + c*n の形式で表現し、係数 a, b, c がどんな値になるかを考察してみた結果だといえる。
考察の際、「組み合わせパターン数」を加味してみた訳である。
振り返って s1(n-1) を同じ形式で表記すると下のようになる。
s1(n-1) = 1/2 * (c(2,0)*B0 * n^2 + c(2,1)*B1 * n)
☆ 注目の式： s3(n-1) = 1/4(n^4 - 6*s2(n-1) - 4*s1(n-1) - n)
s3(n-1) = 1/4(n^4 - 2*n^3 + n^2)
まず、「組み合わせパターン数」(２項係数)をざっとおさらいしておこう。
１乗、２乗、３乗、４乗を記してみる。
(a+b)^1 = a + b
(a+b)^2 = a^2 + 2ab + b^2
(a+b)^3 = a^3 + 3a^2b + 3ab^2 + b^3
(a+b)^4 = a^4 + 4a^3b + 6a^2b^2 + 4ab^3 + b^4
上の展開式の係数だけ取り出すと下のとおり。５乗も追加。
いずれも、上って下りてくるパターンだ。しかも左右対称。これをパスカルの三角形というらしい。
累乗の和に話しを戻そう。
２乗の和 s2(n-1) のときは n^3 に関わる [1, 3, 3, 1] を用いた。
３乗の和 s3(n-1) に対しては [1, 4, 6, 4, 1] を使う。
n^4 = 4 * s3(n-1) + 6 * s2(n-1) + 4 * s1(n-1) + n
上の式が「２乗の和」の式から類推できるだろうか。ピンとこなければ、前の説明を振り返っていただきたい。
ここで s2(n-1) と s1(n-1) は既に分かっている。
しかも、s2(n-1) は s1(n-1) で表現できる。
u = s1(n-1) とした場合、s2(n-1) = 1/3(n^3 - 3u - n) であり、s3(n-1) は次のような式になる。
s3(n-1) = 1/4*(n^4 - 6/3*(n^3 - 3u - n) - 4u - n)
s3(n-1) = 1/4*(n^4 - 2*n^3 + 2u + n)
u = 1/2*(n^2 - n) を当てはめると、結局 下のとおり。
s3(n-1) = 1/4*(n^4 - 2*n^3 + n^2)
「組み合わせパターン数」を持ち込んで式を書き換えると
s3(n-1) = 1/4*(c(4,0)*n^4 + c(4,1)*(-1/2)*n^3 + c(4,2)*(1/6)*n^2)
n^4, n^3, n^2 は出てくるが n が消えてしまっている。
これを無理矢理くっつけるとすれば、c(4,3)*0*n を置くことになるだろうか。
この無理矢理ひねり出した０が、ベルヌーイ数 B3 の値である。
念のため s4(n-1) も見てみよう。
n^5 = 5*s4(n-1) + 10*s3(n-1) + 10*s2(n-1) + 5*s1(n-1) + n
u = s1(n-1) として式を整理すると、5u - 10u + 5u で u が消えてしまい、結局 次のようになる。
s4(n-1) = 1/5*(n^5 - 5/2*n^4 + 5/3*n^3 - 1/6*n)
今度は n^2 が出てこないが、これは B3 = 0 に対応するものである。
s4(n-1) = 1/5*(c(5,0)*n^5 + c(5,1)*(-1/2)*n^4 + c(5,2)*(1/6)*n^3 +
c(5,3)*0*n^2 + c(5,4)*(-1/30)*n)
新たに B4 = -1/30 が判明した。
この先 s5(n-1) を計算してみると分かるが、B5 = 0 となる。
実は B3, B5, B7, B9, …… のいずれも０である。
B1 = -1/2 は例外だが、奇数番号のベルヌーイ数は 0 となる。
これまでのやり方からすると、ベルヌーイ数は、二項定理により得られた展開式からひねり出した「つじつま合わせの数」という印象を受ける。
とはいえ、２乗・３乗・４乗の和のいずれにも規則どおり顔を覗かせる。ただの「いいかげんな数」ではなさそうだ。
それにしても、二項定理を逐次 応用する方法は、入れ子構造が反復して出てくるため面倒だ。
以降では、ベルヌーイ数を算出するための定式化されたノウハウをみていくことにしよう。
ベルヌーイ数を算出する２種類の方法を取り上げる。
B0を土台にしてB1を算出し、次に B0とB1を基にしてB2を算出する随時方式、
B1, B2, B3, B4 など特定の番号のベルヌーイ数を直接割り出す方式の２種類。
それらの算出方法を数式なしで表現するのは難しいので、プログラムで表記する。
また、累乗の和をベルヌーイ数を用いて算出するプログラムも掲げる。
rubyというスクリプト系言語を用いるが、次の２点を念頭に置いておけば容易に理解できるとおもう。
1/2 は Rational(1,2) と表現する。Rational(2) (引数が一つ)は Rational(2,1) とみなされる。2^5 は 2**5 と表現する。
ベルヌーイ数の０番 B0 は、１であることを前提にする。
B1 以降の算出は下のようにして行う。
B1 = -1/2 * c(2,0)*B0 = -1/2
B2 = -1/3 * (c(3,0)*B0 + c(3,1)*B1) = 1/6
B3 = -1/4 * (c(4,0)*B0 + c(4,1)*B1 + c(4,2)*B2) = 0
B4 = -1/5 * (c(5,0)*B0 + c(5,1)*B1 + c(5,2)*B2) = -1/30
…………
Bn = -1/(n+1) * (
c(n+1, 0)*B0 + c(n+1, 1)*B1 + c(n+1, 2)*B2 + ……
c(n+1, (n-1))*B[n-1])
c(n,k) は、n個からk個を取り出す「組み合わせパターン数」を示す。
これをrubyプログラムで書くと下のようになる。
# 組み合わせパターン数(二項係数)を算出
def nck(n, k)
kk = (1..k).to_a # 1, 2, 3, …… k を配列にする
c = Rational(1)
for i in (n-k+1)..n
c *= i
c /= kk.pop
end
return c.to_i
end
# n番目までの一連のベルヌーイ数を算出
def bern1(num)
bn = [] # この配列にベルヌーイ数を入れる
bn[0] = 1
for n in 1..num
s = Rational(0)
for i in 0..(n-1)
s = s + nck(n+1, i) * bn[i]
end
bn[n] = -Rational(1, n+1) * s
end
return bn
end
# テスト表示
bn = bern1(30) # B0〜B30 を取得
for i in 0..30
next if bn[i] == 0 # ベルヌーイ数が0ならスキップ
printf("B%d: %s\n", i, bn[i])
end
中核は def bern1(num) 以降の関数定義部分である。
bn[0] = 1 と定義し、その後は随時 bn[1], bn[2], …… bn[num] を求める。
そこのロジックは、前述のとおり 二項係数とベルヌーイ数を掛け算し、更にそれらを足し算する形で計算を進める。
テスト表示では B30 までを求め、値が０になる B3, B5, B7 などを除外して、ベルヌーイ数を表示している。
rubyの有理数型 Rational は、分子と分母を別々に保持し(既約分数として保持)、どちらもBigNum(巨大な整数)に対応している。
処理時間を気にしないのなら、数値の大きさに関する制約はないと考えていいようだ。
今度は特定の番号のベルヌーイ数を算出する方法である。
B4 を求める手順を見てみよう。
0番目から4番目までの５種類の計算を行い、最後にそれらを合算する。
上の 0番目〜4番目の結果を足し算すると B4 = -1/30 が算出される。
「0番目」も本来は (-1)^0 * 0^4 * …… と計算するのが本筋だが、計算するまでもなく 0 になるので 0 とだけ記した。
ただし、B0 の算出の際は 0^0 が 0 でなく 1 になるので要注意。
n番のベルヌーイ数を算出するプログラムは下のようになる。
# n番目のベルヌーイ数を算出(ベルヌーイ数の一般項)
def bern2(num)
bn = Rational(0) # 最終的にこれにベルヌーイ数が入る
for k in 0..num
elm1 = (-1)**k # (-1)^k
elm2 = k**num # k^num
elm3 = Rational(0) # 以下の総和をこれに代入
for n in k..num
elm3 = elm3 + Rational(1, n+1) * nck(n, k)
end
bn = bn + elm1 * elm2 * elm3
end
return bn
end
# テスト表示
for i in 0..30
bn = bern2(i)
next if bn == 0
printf("B%d: %s\n", i, bn)
end
プログラムとしてはごく簡単だが、手計算で算出するとなると面倒だ。
B30 を求めるときは 30^30 つまり 30の30乗が出てくる。
実際の計算処理はともかくとして、プログラムをみると、計算の手順そのものは単純であることが分かる。
二項定理により累乗の和を計算する際、ベルヌーイ数が顔を覗かせることは既に見たとおりである。
では、ベルヌーイ数を用いて累乗の和を計算するときの手順は……
1^3 + 2^3 + 3^3 + 4^3 + 5^3 = 225 の場合でいうと、
まず下の 0番目から3番目の計算をする。
上の 0番目から3番目までを合算し、それを 3+1 つまり 4 で割り算する。
すると、目的の 225 が算出される。
プログラムとして掲げると次のとおり。
# 1^k + 2^k + …… + n^k を算出：ベルヌーイ数を利用
def bsum(n, k)
bn = bern1(k)
n += 1
s = Rational(0)
for i in 0..k
s = s + nck(k+1, i) * bn[i] * n**(k-i+1)
end
s = s / (k+1)
return s.to_i
end
s = bsum(5, 3) とすれば 5^3 までの和 225 が変数 s に代入される。
zip圧縮ファイルに入っているサンプルプログラムでは、単純な累乗の和の計算を行う nsum() という関数も設け、bsum() と値が一致するか確認できるようになっている。
ベルヌーイ数は、無限級数とも関連が深いようだ。
なぜ級数にベルヌーイ数が登場するのか、私には納得するための道筋が分からないままだが、どんな形かを少しだけ見てみたい。
x / (e^x - 1) の展開式にベルヌーイ数が現れることは多くの解説に出てくる。
というより、x / (e^x - 1) の係数でベルヌーイ数を定義するらしい。
x / (e^x - 1) =
B0/0! * x^0 +
B1/1! * x^1 +
B2/2! * x^2 +
B3/3! * x^3 +
B4/4! * x^4 + ……
上の式が成り立つとのことである。ほんとだろうか。
まず e^x の展開式だが、とても調和のとれた形をしている。
e^x = 1/0! + 1/1!x + 1/2!x^2 + 1/3!x^3 + 1/4!x^4 + ……
何度 微分しても変化しない展開式である。
’e^x - 1’ は、この展開式から定数項の１を引き算したものだ。
この定数項なしの展開式について、その逆数に当たる式を求めるのは難しい。
方程式 f(x) の逆数に相当する 1/f(x) を求めるには、
求める式を a0 + a1x + a2x^2 + a3x^3 + …… と仮定し、
f(x) と掛け算した結果が 1 になるように a0, a1, a2, a3, …… を決めればいい。
あるいは、1/(1-r) = 1 + r + r^2 + r^3 + …… を利用し、
r = 1 - f(x) と置いて計算すれば 1/f(x) が求まる。
実際に計算してみれば分かるが、どちらの方法を採るにしても f(x) の定数項が 1 でないと計算に難渋する。
というより、そもそも計算できないか、あるいは、意味のある結果が得られないとおもう。
e^x - 1 に話を戻そう。
e^x - 1 は定数項なしになってしまったが、
これを x で割り算すれば再び定数項が 1 になる。
しかも、x で割り算した結果の逆数をとれば、
それが求めたい x / (e^x - 1) になる。
こうして求めた展開式に対して、「定数項を記録しては微分する」という処理を繰り返せば、その記録がベルヌーイ数になっているはずだ。
この手順で計算することにしよう。
ここでは、パソコンで処理することを考える。
とはいえ、何もないところからプログラムを作るのは大変なので、拙作 binomer.rb を利用する。
binomerについては ２項定理とそのプログラミング処理・追加編 を参照されたい。
ともあれプログラムを掲げてみる。
require "./binomer" # binomer.rb を取り込む
max = 16 # 無限級数を何項まで扱うかの数
es = Binomer::e(max) # e^x の展開式
w = (es - 1) / "x" # (e^x - 1) / x
bs = w.reciprocal(max) # 逆数をとり x / (e^x - 1) を得る
printf("'x/(e^x-1)':\n%s\n\n", bs.cx.join(" + ")) # x/(e^x-1) を表示
bn = [] # この配列にベルヌーイ数を入れる
while a = bs.first # 初項の情報を取得
if a[2] >= 1 then # 定数項がない(x^0がなく xの指数が1以上)
bn << 0 # ベルヌーイ数は０
else # 定数項があるので
bn << a[0] # 定数項をベルヌーイ数として記録
end
bs = bs.differential # 級数を微分する
end
puts "'Bernoulli number':"
for i in 0...bn.size
next if bn[i] == 0
printf("B%d: %s\n", i, bn[i])
end
binomerでは初等関数 e, log, sin, cos の級数が予め用意されている。
es = Binomer::e(16) とすれば、変数 es に e^x の展開式が第16項まで記録される。
es には 16項の情報が含まれる訳だが、一つの「項」について
係数およびxの指数が数値として保持されている。
a = es.first とすれば、変数 a に展開式の初項の情報が入る。
変数 a は配列で、a[0] に係数、a[2] にxの指数が代入されている。
a[2] が0なら定数項であり、1以上だと x^1 とか x^2 などの「項」である。
なお、a[1] は通常 1 なので無視してかまわない。
次の３行で、目的の x / (e^x - 1) の展開式を変数 bs に代入できる。
es = Binomer::e(16) # e^x の展開式
w = (es - 1) / "x" # (e^x - 1) / x
bs = w.reciprocal(16) # 逆数をとり x / (e^x - 1) を得る
あとは bs の初項から情報を取り出し、微分したら再び初項を見て …… という処理を繰り返す。
繰り返しは while ループで行っている。
なお、第何項まで取り扱うかの数 16 を大きくする場合は注意していただきたい。
Windows7(64bit) の私の環境では、16だと 2秒弱、20の場合は 8秒程度、30にすると 3分11秒かかった。
binomerでは計算量を抑制するための工夫をしていないため、そうした結果になってしまうのだとおもう。
本論とは関係ないプログラムの話が長くなってしまった。この辺にしよう。
ζ(2) が π^2/6 であることを示したのはオイラーだが、彼は ζ(4), ζ(6) などの引数が偶数の場合についていくつか計算したらしい。
実は、ζ(偶数) の値にはベルヌーイ数が関係している。
ζ(2n) = 2^(2n-1) * π^(2n) * |B[2n]| / (2n)!
|B[2n]| は、2n番目のベルヌーイ数の絶対値を示す。
具体的には次のような値になる。
ζ(2) = π^2/6
ζ(4) = π^4/90
ζ(6) = π^6/945
オイラーが上の ζ(2n) の式を導き出していたのかどうかは分からない。
ここで、リーマンの ζ(s) と ζ(1-s) の関係式の適用を考えてみよう。
関係式を使えば、ζ(2) から ζ(-1) が、ζ(4) から ζ(-3) が示せるはずだ。
だとすれば、ζ(マイナス奇数) がベルヌーイ数に関連することになる。
リーマンの関係式は次のとおり。
ζ(1-s) = ζ(s) * 2 * (2π)^(-s) * Γ(s) * cos(s*π/2)
Γ(s) は、s が自然数の場合 (s-1)! つまり (s-1) の階乗になる。
ζ(s) のところに先の ζ(2n) を当てはめると、πについては π^s と π^(-s) が出てくるため、打ち消しあってなくなってしまう。
整理すると下のような式になる。
ζ(1-2n) = |B[2n]| / (2n) * cos(n*π)
ここで cos(n*π) は、nが奇数なら -1、偶数だと 1 になる。プラスとマイナスの記号が切りかわるスイッチの役割といっていい。
なので、ベルヌーイ数に付いている絶対値の記号をはずし、 -1 を掛け算する形にすれば、cos() の代替になる。
結局、下のような簡単な式になる。
ζ(1-2n) = B[2n] / (2n) * -1
具体例をあげると次のとおり。
ζ(-1) = (1/6) / 2 * -1 = -1/12
ζ(-3) = (-1/30) / 4 * -1 = 1/120
ζ(-5) = (1/42) / 6 * -1 = -1/252
なぜベルヌーイ数がこんな舞台に登場するのか、よく分からない。
というより、自然対数の底eが絡む式にしても、リーマンの関係式にしても、なぜそのような式が出てくるのか、私にはその訳が分かっていない。まだまだ精進が足りない。
ベルヌーイ数については他にもいろんな話題がありそうだが、とりあえずの初歩ということで、この辺でおわりにしよう。
〜 「ベルヌーイ数が顔をのぞかせるとき」おわり 〜
Copyright (C) T. Yoshiizumi, 2015 All rights reserved.