Essays/Bernoulli Numbers

From J Wiki
Jump to: navigation, search

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 1

GKP 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*n

Since 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 Roger Hui.