書いていたプログラムでベルヌーイ数が必要になった。
ただ、必要な範囲が決まっているのでわざわざその都度ベルヌーイ数を求めるよりもいっその事最初から配列に格納しておいてそれを使う方が楽だし速い。メモリが限られてるわけでもないし。
ただ、某百科事典には29番目までしか載ってなく、それよりも大きな値も欲しかった。
そこで、ベルヌーイ数を求めるプログラムを書いてみた。
ベルヌーイ数の求め方
以下の漸化式で簡単に求められます。
\begin{align}
B_0 &= 1 \\
B_n &= -\frac{1}{n + 1}\sum_{n = 0}^\infty\hspace{-10pt}\quad_{n + 1}\mathrm{C}_kB_k
\end{align}
プログラム
Python3
Python 3.6.4で動きました。なお、僕は普段Pythonを使わないのですが、
"PythonならInt型でオーバーフローしないじゃん"
と思ったので慣れないPythonを使いました。
(ちなみに僕は動的型付け言語が大嫌いなので必要にならない限り使わない)
ソースコード
def gcd(arg1, arg2):
if arg2 > arg1:
return gcd(arg2, arg1)
while True:
tmp = arg2
arg2 = arg1 % arg2
arg1 = tmp
if arg2 == 0:
break
return arg1
def lcm(arg1, arg2):
return (arg1 * arg2) // gcd(arg1, arg2)
def r_add(rval1, rval2):
s = lcm(rval1[1], rval2[1])
return (rval1[0] * (s // rval1[1]) + rval2[0] * (s // rval2[1]), s)
def r_mul(rval1, rval2):
return (rval1[0] * rval2[0], rval1[1] * rval2[1])
def r_reduce(rval):
if rval[0] == 0:
return rval
s = gcd(rval[0], rval[1])
r = (rval[0] // s, rval[1] // s)
if r[1] < 0:
r = (-r[0], -r[1])
return r
def combination(arg1, arg2):
if arg1 - arg2 < arg2:
arg2 = arg1 - arg2
r = 1
for i in range(arg2):
r *= arg1 - i
for i in range(arg2):
r //= i + 1
return r
n = 60
B = [(1, 1), (-1, 2)]
for i in range(2, n + 1):
if i & 1 != 0:
B.append((0, 1))
continue
s = (0, 1)
for j in range(i):
if B[j][0] == 0:
continue
c = combination(i + 1, j)
s = r_add(s, r_mul((c, 1), B[j]))
r = r_mul(s, (-1, i + 1))
B.append(r_reduce(r))
for i in range(0, n + 1):
if B[i][0] != 0:
print("B{0} = {1} / {2} := {3:.2e}".format(i, B[i][0], B[i][1], B[i][0] / B[i][1]))
else:
print("B{0} = 0".format(i))
分数はタプルで表しており、$\frac{1}{2}$は(1, 2)と表しています。
そして、タプルで表された分数の和をとるr_add関数、積をとるr_mul関数、約分をするr_reduce関数があり、それらの関数で使うために整数の最大公約数と最小公倍数を求めるgcd,lcm関数があります。
また、漸化式で使うのでcombination関数も定義してあります。
あとは漸化式を計算してリストBに格納し、最後に順番に表示しているだけです。
何番目まで計算するかは変数nによって変わります。n = 60となっているのでこれは$B_0$から$B_{60}$まで求めます。適宜nの値変えてください。
ちなみに、$B_0$は漸化式の初期値なので入れるのは当たり前として、$B_1$の値もプログラムに直接書き込んでいます。ベルヌーイ数は$B_1$以外の奇数番目の値が0になるので、こうする事で奇数番目の処理が楽になるからです。
結果
B0 = 1 / 1 := 1.00e+00
B1 = -1 / 2 := -5.00e-01
B2 = 1 / 6 := 1.67e-01
B3 = 0
B4 = -1 / 30 := -3.33e-02
B5 = 0
B6 = 1 / 42 := 2.38e-02
(以下略)
ベルヌーイ数 一覧表
$n$ | 分子 | 分母 | 実数値 |
---|---|---|---|
0 | 1 | 1 | 1.00e+00 |
1 | -1 | 2 | -5.00e-01 |
2 | 1 | 6 | 1.67e-01 |
3 | 0 | - | 0 |
4 | -1 | 30 | -3.33e-02 |
5 | 0 | - | 0 |
6 | 1 | 42 | 2.38e-02 |
7 | 0 | - | 0 |
8 | -1 | 30 | -3.33e-02 |
9 | 0 | - | 0 |
10 | 5 | 66 | 7.58e-02 |
11 | 0 | - | 0 |
12 | -691 | 2730 | -2.53e-01 |
13 | 0 | - | 0 |
14 | 7 | 6 | 1.17e+00 |
15 | 0 | - | 0 |
16 | -3617 | 510 | -7.09e+00 |
17 | 0 | - | 0 |
18 | 43867 | 798 | 5.50e+01 |
19 | 0 | - | 0 |
20 | -174611 | 330 | -5.29e+02 |
21 | 0 | - | 0 |
22 | 854513 | 138 | 6.19e+03 |
23 | 0 | - | 0 |
24 | -236364091 | 2730 | -8.66e+04 |
25 | 0 | - | 0 |
26 | 8553103 | 6 | 1.43e+06 |
27 | 0 | - | 0 |
28 | -23749461029 | 870 | -2.73e+07 |
29 | 0 | - | 0 |
30 | 8615841276005 | 14322 | 6.02e+08 |
31 | 0 | - | 0 |
32 | -7709321041217 | 510 | -1.51e+10 |
33 | 0 | - | 0 |
34 | 2577687858367 | 6 | 4.30e+11 |
35 | 0 | - | 0 |
36 | -26315271553053477373 | 1919190 | -1.37e+13 |
37 | 0 | - | 0 |
38 | 2929993913841559 | 6 | 4.88e+14 |
39 | 0 | - | 0 |
40 | -261082718496449122051 | 13530 | -1.93e+16 |
41 | 0 | - | 0 |
42 | 1520097643918070802691 | 1806 | 8.42e+17 |
43 | 0 | - | 0 |
44 | -27833269579301024235023 | 690 | -4.03e+19 |
45 | 0 | - | 0 |
46 | 596451111593912163277961 | 282 | 2.12e+21 |
47 | 0 | - | 0 |
48 | -5609403368997817686249127547 | 46410 | -1.21e+23 |
49 | 0 | - | 0 |
50 | 495057205241079648212477525 | 66 | 7.50e+24 |
51 | 0 | - | 0 |
52 | -801165718135489957347924991853 | 1590 | -5.04e+26 |
53 | 0 | - | 0 |
54 | 29149963634884862421418123812691 | 798 | 3.65e+28 |
55 | 0 | - | 0 |
56 | -2479392929313226753685415739663229 | 870 | -2.85e+30 |
57 | 0 | - | 0 |
58 | 84483613348880041862046775994036021 | 354 | 2.39e+32 |
59 | 0 | - | 0 |
60 | -1215233140483755572040304994079820246041491 | 56786730 | -2.14e+34 |
雑記
僕はWordが嫌いなのでLatexをよく使います。なので、Qiitaの投稿でLatexの数式が使えるのはすごく便利なのですが、いまいちカバーしきれてないのが難点です。
漸化式を書くのに組み合わせCの左下の添字を書こうとして
\sum_{n = 0}^\infty{}_{n + 1}\mathrm{C}_kB_k
と書いたら数式が表示されなかった。
ちなみにLaTeXiTでの表示結果は
こうなるのであってます。
でも、左下の添字が使えず、なぜか左上のみの添字と、左上と左下両方の添字なら使える。
なので、苦肉の策として
\sum_{n = 0}^\infty\hspace{-10pt}\quad_{n + 1}\mathrm{C}_kB_k
左上に空白の添字を追加して、そうすると不自然に空白ができるのでマイナスにした\hspaceで調整するということをした。
でも、この使いづらさ、なんとかならないのかな...