ベルヌーイ数の計算
calculation of Bernoulli numbers
世界初のコンピュータープログラムは、ベルヌーイ数を求めるものであったとされている。
そのベルヌーイ数を分数として、漸化式を使用せずに直接求めるプログラムを考える。
二項係数を使用しないベルヌーイ数の一般式は (ベルヌーイ数の生成規則(1) を参照)
( k=1,2,3,4,… )
ベルヌーイ数を四則とべき乗演算を使用して求める関数のUBASICプログラム。奇数の場合も計算値を返している。
800 ' function : Bernoulli number
(0<= K <=520)
810 fnB(K)
820 local Bn=0^K,F=1,I,J,M,R,S
830 for J=1 to K:S=0
840 for M=1 to K-J+1
850 R=1:for I=M to
M+J-1:R*=I:next:S+=R//I
860 next
870 F*=-J:Bn+=J^K//F*S
880 next
890 return(Bn)
?.B(0),.B(1),.B(2),.B(3),.B(4),.B(5),.B(6),.B(8),.B(12),.B(80)
1 -1//2
1//6 0 -1//30
0 1//42 -1//30 -691//2730
-4603784299479457646935574969019046849794257872751288919656867//230010
このプログラムは処理速度は遅いが、全て分数計算のため値は正確である。ベルヌーイ数そのものを求めるのに適している。実際に求めた値は次の通り。
0
1 /
1
1.0000000000…
2
1 /
6
0.1666666666…
4
-1 / 30 -0.0333333333…
6
1 /
42
0.0238095238…
8
-1 /
30
-0.0333333333…
10
5 /
66
0.0757575757…
12
-691 /
2730
-0.2531135531…
14
7 /
6
1.1666666666…
16
-3617 /
510
-7.0921568627…
18
43867 /
798
54.9711779448…
20
-174611 /
330
-529.1242424242…
22
854513 /
138
6192.1231884057…
24
-236364091 /
2730
-86580.2531135531…
26
8553103 /
6
1425517.1666666666…
28
-23749461029 /
870
-27298231.0678160919…
30
8615841276005 /
14322
601580873.9006423683…
32
-7709321041217 /
510
-15116315767.0921568627…
34
2577687858367 /
6
429614643061.1666666666…
36 -26315271553053477373 /
1919190 -13711655205088.3327721590…
38 2929993913841559 /
6
488332318973593.1666666666…
40 -261082718496449122051 /
13530 -19296579341940068.1486326681…
二項係数を使用する次のベルヌーイ数の一般式によるプログラム。
( k=0,1,2,3,… )
800 ' function : Bernoulli number
(0<= K <=699)
810 fnB(K)
820 local Bn,D,J,M
830 for J=0 to K
840 D=0:for M=J to K:D+=combi(M,J)//(M+1):next
850 Bn+=(-1)^J*D*J^K
860 next
890 return(Bn)
?.B(0),.B(1),.B(2),.B(3),.B(4),.B(5),.B(6),.B(8),.B(12),.B(80)
1 -1//2
1//6 0 -1//30
0 1//42 -1//30 -691//2730
-4603784299479457646935574969019046849794257872751288919656867//230010
最後以外は整数計算である次のベルヌーイ数の一般式によるプログラム。
( k=1,2,3,4,… )
タンジェント数 |T| を同時に求めることができる。
800 ' function : Bernoulli number
(0<= k <=853)
810 fnB(K)
820 local Bn=1,D,J,M,T
830 for M=1 to K
840 D=0:for J=M to
K:D+=combi(K,J):next
850
T+=(-1)^(M-1)*D*M^(K-1)
860 next
870 if K>0 then
Bn=(-1)^K*K*T//(2^K*(2^K-1))
890 return(Bn)
?.B(0),.B(1),.B(2),.B(3),.B(4),.B(5),.B(6),.B(8),.B(12),.B(80)
1
-1//2 1//6 0
-1//30 0 1//42
-1//30 -691//2730
-4603784299479457646935574969019046849794257872751288919656867//230010
項数の少ないベルヌーイ数の一般式による処理が速いプログラム。
(
k=2,3,4,5,… ) 
( k=0,1
)
重複計算を避けるために次の例のように順序を変更し、二項係数を乗算で計算する。
B4= {(-1/2)4} { (1)33- (1+(4/1))23 + (1+(4/1)+(4/1)(3/2))13 } {4/(24-1)}
800 ' function : Bernoulli number
(0<= k <=853)
810 fnB(K)
820 local Bn=(-1//2)^K,C=1,D=1,J,M,T
830 for M=K-1 to 1 step -1
840 T+=(-1)^(M-1)*D*M^(K-1)
850 C*=(M+1)//(K-M):D+=C
860 next
870 if K>1 then Bn*=T*K//(2^K-1)
890 return(Bn)
?.B(0),.B(1),.B(2),.B(3),.B(4),.B(5),.B(6),.B(8),.B(12),.B(80)
1 -1//2
1//6 0 -1//30
0 1//42 -1//30 -691//2730
-4603784299479457646935574969019046849794257872751288919656867//230010
範囲が広くて以上の中では最速のプログラム。奇数のタンジェント数の式を利用する。
( k=0,1,2,3,… )
T5= (1)35 - (1+(7/1))25 + (1+(7/1)+(7/1)(6/2))15
B6= (-1)6/2-16T5/(26(26-1)) のように計算する。実際の計算方法は少し異なる。
800 ' function : Bernoulli
number (0<= k <=955)
810 fnB(K)
820 local
Bn=1,C=1,D=1,J,M,N=K//2,T
830 if K=1 then
return(-1//2) else if odd(K) then return(0)
840 for M=N to 1
step -1
850
T+=(-1)^M*D*M^(K-1)
860
C*=(N+1+M)//(N+1-M):D+=C
870 next
880 if K>1
then Bn=-T*K//(2^K*(2^K-1))
890 return(Bn)
?.B(0),.B(1),.B(2),.B(3),.B(4),.B(5),.B(6),.B(8),.B(12),.B(80)
1 -1//2
1//6 0 -1//30
0 1//42 -1//30 -691//2730
-4603784299479457646935574969019046849794257872751288919656867//230010
クラウゼン・フォンシュタウトの定理(Clausen-von Staudt's theorem)によると、ベルヌーイ数は素数の逆数の和と整数に分解できる。
Bk=Ck-Σp-1|k 1/p ( kは1と2,4,6,… Ckは整数,pは素数, kはp-1で割り切れる )
2以上の偶数の場合を次表に示す。
|
B2k |
値 |
整数-Σ(1/素数) |
B2k |
値 |
整数-Σ(1/素数) |
|
B2 |
1/6 |
1-1/2-1/3 |
B12 |
-691/2730 |
1-1/2-1/3-1/5-1/7-1/13 |
|
B4 |
-1/30 |
1-1/2-1/3-1/5 |
B14 |
7/6 |
2-1/2-1/3 |
|
B6 |
1/42 |
1-1/2-1/3-1/7 |
B16 |
-3617/510 |
-6-1/2-1/3-1/5-1/17 |
|
B8 |
-1/30 |
1-1/2-1/3-1/5 |
B18 |
43867/798 |
56-1/2-1/3-1/7-1/19 |
|
B10 |
5/66 |
1-1/2-1/3-1/11 |
B20 |
-174611/330 |
-528-1/2-1/3-1/5-1/11 |
またkが偶数の場合 Bk=-2ζ(k)k!/(2πi)k で、Ckはこの値に近い整数である。k<720 では
となる(下記参照)。次のプログラムではζ(k)を少ない回数の計算で置き換えて、ベルヌーイ数を求めている。
(k=2,4,6,…)
900 ' function : Bernoulli number
910 fnB(K)
920 local A,Bn,D,J,P=2
930 if K=1 then return(-1//2) else
if odd(K) then return(0)
940 while (P<K+1)
950 inc
J:P=prm(J):D+=1/J^K:A+=(K@(P-1)=0)//P
960 wend
980
Bn=int(pi(1/2)-D*2*!(K)/pi(2#i)^K)-A
990 return(Bn)
?.B(0),.B(1),.B(2),.B(3),.B(4),.B(5),.B(6),.B(8),.B(12),.B(80)
1
-1//2 1//6 0
-1//30 0 1//42
-1//30 -691//2730
-4603784299479457646935574969019046849794257872751288919656867//230010
odd() は奇数かどうかの関数、prm(J) はJ番目の素数、@ は余り演算、pi(c)はcπ を表す。
(K@(P-1)=0)の値は余りが0ならば1(真)で、余りが0以外ならば0(偽)となる。
このプログラムは繰り返しが(k以下の素数の個数)+1回以下であるので処理が速い。
ただしπk等を使用した無理数の計算であるため、小数点以下の桁数(UBASICではpoint文で指定する)により精度が決まる。例えばpoint 15 は10進数で72桁の指定であるが、この場合はB92までは正確な値を返す。またpoint 150 ではB492まで正確である。この近辺以上では整数部分の容量不足によりオーバーフローが発生するため、最大値が低下する。
無理数等が混在する式で、ベルヌーイ数を過剰な精度ではなく、実数と同じレベルの精度で速く求めて使用する目的に適している。
point 小数
正確 point 小数 正確 point 小数
正確 point 小数 正確
15 72
92 55 264
230 95 457
346 135 650 452
20
96 112 60 288
246 100 481
360 140 674 466
25 120
132 65 313
260 105 505
374 145 698 478
30
144 148 70
337 274 110
529 388 150
722 492
35 168
166 75 361
290 115 553
400 155 746 490
40
192 182 80
385 304 120
577 414 160
770 472
45 216
198 85 409
318 125 602
428 165 794 454
50
240 214 90 433
332 130 626
440 170 818 434
誤差の程度としてB94の例を示す。
上側がpoint 150 で計算した正確な値で、下側がpoint 15 で計算した値である。
1220813806579744469607301679413201203958508415202696621436215105284649447/6
1220813806579744469607301679413201203958508415202696621436215105284649459/6
分数を実数に変換した値は
203468967763290744934550279902200200659751402533782770239369184214108241.16666…
203468967763290744934550279902200200659751402533782770239369184214108243.16666…
相対誤差は2/B94であり、point 15(小数点以下72桁)の表示では
0.000000000000000000000000000000000000000000000000000000000000000000000009
偶数ゼータ関数値を求めるプログラムは次の通り。
900 ' function : zeta(2n)
910 fnZeta(K)
920 local A,Bn,D,F,J,P=2,Zt
930 if odd(K) then stop:else if K<0 then return(0)
940 while (P<K+1)
950 inc J:P=prm(J):D+=1/J^K:A+=(K@(P-1)=0)//P
960 wend
970 F=2*!(K)/pi(2#i)^K
980 Bn=int(pi(1/2)-D*F)-A:Zt=-Bn/F
990 return(Zt)
ζ(2)/22+ζ(4)/24+ζ(6)/26+…=1/2 の式で確認する。
point -50:a=0:for n=1 to
400:a+=.zeta(2*n)/2^(2*n):next:?a
0.49999999999999999999999999999999999999999999999999999999999999999999999999999
99999999999999999999999999999999999999999999999999999999999999999999999999999999
99999999999999999999999999999999999999999999999999999999999999999999999999999999
989
更に …+ζ(-4)/2-4+ζ(-2)/2-2+ζ(0)/20+ζ(2)/22+ζ(4)/24+…=0 の式で確認する。
point -50:a=0:for n=-400 to 400:a+=.zeta(2*n)/2^(2*n):next:?a
-0.00000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000000000000000
01
Bk=Ck-Σp-1|k 1/p を Bk=Ck-Σ1/p と略記する。Σ1/p はkが大きくなるにつれて徐々に最大値が大きくなるが、速度は非常に遅い。そのため整数 Ck=Bk+Σ1/p の値を
で求めることが可能な条件を調べる。これがCkと等しいためには 0<=f-Σ1/p<1 つまり Σ1/p<=f<1+Σ1/p となる必要がある。
Σ1/p の最小値は 1/2+1/3=5/6=0.83333… である。Σ1/p が新しく最大となるkの値を示す。
k Σ1/p
k Σ1/p
k Σ1/p
2 0.83333
1440 1.58923 55440 1.77358
4 1.03333
2160 1.59640 110880 1.78965
12 1.25311
2520 1.60830 166320 1.79644
36 1.33277
5040 1.68408 221760 1.79947
60 1.39267 10080
1.69667 277200 1.79997
120 1.41706 15120
1.70572 332640 1.81258
180 1.47786 20160
1.70586 554400 1.81605
240 1.48004 25200
1.70936 665280 1.82266
360 1.51595 30240
1.71834 720720
1.83497
720 1.57892
50400 1.72195
: :
従って 0.83333…<=f<1.83333… となる必要があるが、なるべく広い範囲で成立させるためにはfの値を大きくしなければならない。
上記式の f=(1/1-1/2)π=π/2=1.57079… の場合はkが720未満で成立する。
f=(1/1-1/2+1/3-1/4)π=7π/12=1.83259… の場合はkが720720未満で成立し、範囲が約1,000倍広くなる。
ベルヌーイ数を階乗と自然対数の値により求める。
級数と乗積の「自然対数と階乗・整数」と「自然数の累乗の階差とゼータ関数」を参照。
210
point -150:F=10^18
220 for N=0 to
17
230 U=round(!(N)/log(1+1/F)^(N+1)):B=0
240
for K=1 to N+1
250 W=U/(-F)^(K+1)
260 B+=round((W-round(W))*F)//K
270
next
280
print (-1)^N*B,
290 next
1
-1//2 1//6 0
-1//30 0 1//42
0 -1//30 0
5//66 0
-691//2730
0 7//6
0 -3617//510
0