Overview
-The Bernoulli numbers could be computed by the relations:
B(0) = 1 ; B(0) + Cn+11 B(1) + Cn+12 B(2) + ...... + Cn+1n B(n) = 0 where Cnk = n!/(k!(n-k)!) are the binomial coefficients
-However, this recurrence relation is unstable and the results are quite
inaccurate for large n ( for n = 20 , only 4 digits are correct!
)
-The following program uses a series expansion instead:
B(n) = (-1)-1+n/2 2n!/(2pi)n ( 1/1n + 1/2n + ...... + 1/kn + ...... ) if n is even and B(0) = 1 ; B(1) = -1/2 ; B(2n+1) = 0 if n > 0
-Actually, B(2) = 1/6 ; B(4) = -1/30 ; B(6) = -1/42 are given
directly ( lines 32 to 39 may be deleted without a great
loss of speed )
Program listing
Data Registers: R00 to R02: temp
Flags: /
Subroutine: "ZETA" ( cf "Zeta Function
for the HP-41" )
| 01 LBL "BERN"
02 STO 02 03 1 04 X>Y? 05 RTN 06 ST+ X 07 X<=Y? 08 GTO 00 09 1/X 10 CHS 11 RTN 12 LBL 00 13 X#Y? |
14 GTO 00
15 6 16 1/X 17 RTN 18 LBL 00 19 MOD 20 0 21 X#Y? 22 RTN 23 4 24 RCL 02 25 X#Y? 26 GTO 00 |
27 30
28 1/X 29 CHS 30 RTN 31 LBL 00 32 6 33 X#Y? 34 GTO 00 35 42 36 1/X 37 RTN 38 LBL 00 39 X<>Y |
40 FIX 9
41 XEQ "ZETA" 42 ST+ X 43 1 44 CHS 45 RCL 02 46 2 47 / 48 Y^X 49 * 50 CHS 51 PI 52 ST+ X |
53 E-9
54 - 55 RCL 02 56 Y^X 57 / 58 LBL 01 59 RCL 02 60 * 61 DSE 02 62 GTO 01 63 END |
( 91 bytes / SIZE 003 )
| STACK | INPUTS | OUTPUTS |
| X | n | B(n) |
Example: 116 XEQ "BERN"
gives B(116) = -1.748892190 1098
( in 24 seconds )
References:
[1] John H. Conway & Richard K. Guy , "The Book of Numbers"
- Springer Verlag - ISBN 0-387-97993-X
[2] Abramowitz and Stegun , "Handbook of Mathematical Functions"
- Dover Publications - ISBN 0-486-61272-4