LoginSignup
2
0

More than 5 years have passed since last update.

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

Last updated at Posted at 2018-03-09

書いていたプログラムでベルヌーイ数が必要になった。
ただ、必要な範囲が決まっているのでわざわざその都度ベルヌーイ数を求めるよりもいっその事最初から配列に格納しておいてそれを使う方が楽だし速い。メモリが限られてるわけでもないし。
ただ、某百科事典には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での表示結果は
latex-image-1.png
こうなるのであってます。
でも、左下の添字が使えず、なぜか左上のみの添字と、左上と左下両方の添字なら使える。
なので、苦肉の策として
\sum_{n = 0}^\infty\hspace{-10pt}\quad_{n + 1}\mathrm{C}_kB_k
左上に空白の添字を追加して、そうすると不自然に空白ができるのでマイナスにした\hspaceで調整するということをした。
でも、この使いづらさ、なんとかならないのかな...

リンク

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

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0