ベルヌーイ数を分数で求める

書いていたプログラムでベルヌーイ数が必要になった。
ただ、必要な範囲が決まっているのでわざわざその都度ベルヌーイ数を求めるよりもいっその事最初から配列に格納しておいてそれを使う方が楽だし速い。メモリが限られてるわけでもないし。
ただ、某百科事典には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

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

Scroll to top