Bernoulli numbers are defined by the recurrence relation
(0=m) = +/ (B 1+m) * (i.1+m)!1+m
for m>:0 (Graham, Knuth, and Patashnik, Concrete Mathematics, Section 6.5). They can be computed as follows:
B0=: 3 : 0 b=. ,1x for_m. }.i.x: y do. b=. b,(1+m)%~-+/b*(i.m)!1+m end. )
B0 y are the first y Bernoulli numbers. For example:
B0 1 1 B0 13 1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0 5r66 0 _691r2730 (B0 13) +/ .* (i.!]) 13x 0
The number of iterations can be halved by exploiting the fact that the odd-numbered Bernoulli numbers (other than the first) are 0.
B=: 3 : 0 b=. 1 _1r2 for_m. 2x*}.i.>.-:y do. b=. b,0,~(1+m)%~-+/b*(i.m)!1+m end. }:^:(2|y) b ) B 13 1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0 5r66 0 _691r2730 (B +/ .* i.!])"0 >:i.50x 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
A shorter (but less efficient) expression derives by considering the recurrence relation as a matrix equation.
(i.!])"0 >: i.10x
1 0 0 0 0 0 0 0 0 0
1 2 0 0 0 0 0 0 0 0
1 3 3 0 0 0 0 0 0 0
1 4 6 4 0 0 0 0 0 0
1 5 10 10 5 0 0 0 0 0
1 6 15 20 15 6 0 0 0 0
1 7 21 35 35 21 7 0 0 0
1 8 28 56 70 56 28 8 0 0
1 9 36 84 126 126 84 36 9 0
1 10 45 120 210 252 210 120 45 10
(10{.1) %. (i.!])"0 >: i.10x
1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0
B 10
1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0
B1=: {.&1 %. (i.!])@>:@i.@x:
B1 10
1 _1r2 1r6 0 _1r30 0 1r42 0 _1r30 0
(B -: B1)"0 >: i.50
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1GKP equation 6.89 relates the 2*n-th Bernoulli number and the generalized harmonic number H2n of order 2*n :
H2n = (_1^n-1) * (2^(2*n)-1) * (1p1^2*n) * ({:B 1+2*n) % !2*nSince H2n tends to 1 , we can derived from this approximations for the sizes of the Bernoulli numbers and for the log of their absolution values. Thus:
BH=: 3 : 0
if. 2|y do. 0 else. (_1^<:-:y)*+:(!y)%2p1^y end.
)
BL=: 3 : 0
if. 2|y do. __ else. (^.2)+(+/^.>:i.y)-y*^.2p1 end.
)
0j_16 ": BH 100
_2.8382249570693794e78
0j_16 ": {: B 101
_2.8382249570693707e78
^. {: B 171
394.827
BL 170
394.827
Collected Definitions
B0=: 3 : 0
b=. ,1x
for_m. }.i.x: y do. b=. b,(1+m)%~-+/b*(i.m)!1+m end.
)
B=: 3 : 0
b=. 1 _1r2
for_m. 2x*}.i.>.-:y do. b=. b,0,~(1+m)%~-+/b*(i.m)!1+m end.
}:^:(2|y) b
)
B1=: {.&1 %. (i.!])@>:@i.@x:
BH=: 3 : 0
if. 2|y do. 0 else. (_1^<:-:y)*+:(!y)%2p1^y end.
)
BL=: 3 : 0
if. 2|y do. __ else. (^.2)+(+/^.>:i.y)-y*^.2p1 end.
)
Contributed by RogerHui.