Bernoulli Numbers

[Namir]

Bernoulli Numbers using Series Approximations

Algorithm

Given integer n and tolerance toler

Code:

IF n=1 THEN Return -0.5
IF n Mod 2 = 1 THEN Return 0 

sum=1 
i=2 
term=2*toler 
WHILE term>toler DO
  term=1/i^n 
  sum=sum+term 
  i=i+1 
END 
IF n MOD 4 == 0 THEN chs=-1 ELSE chs=1
t=chs*(n!)/pi^n/2^(n-1) 
Return t*sum

HP-41C Implementation

Memory Map


R00 = n
R01 = toler
R02 = i
R03 = Sum
R04 = term


Implementation

Code:

1 LBL "BERSER"
2 LBL A
3 CF 22
4 "TOLER?"
5 PROMPT
6 FS? 22
7 STO 01
8 "N?"
9 PROMPT
10 STO 00
11 1
12 X<>Y
13 X=Y?
14 GTO 11
15 2
16 MOD
17 1
18 X=Y?
19 GTO 10
20 1
21 STO 03      # sum = 1
22 STO 02      # i = 1
23 LBL 01      # Start summation loop
24 1
25 STO+ 02     # i=i+1
26 RCL 02
27 RCL 00
28 Y^X       
29 1/X         # Calculate 1/i^n
30 STO+ 03     # sum = sum + 1/i^n
31 RCL 01
32 X<=Y?       # toler <= term?
33 GTO 01      # resume to the next iteration
34 CF 00       # start calculating other factor
35 RCL 00
36 4
37 MOD
38 X=0?
39 SF 00
40 RCL 00
41 FACT        # n!
42 PI
43 RCL 00
44 Y^X         # pi^n
45 /
46 2
47 RCL 00
48 1
49 -
50 Y^X          # 2^(n-1)
51 /
52 RCL 03
53 *
54 FS?C 00     # change sign?
55 CHS
56 RTN
57 LBL 10
58 0
59 RTN
60 LBL 11
61 -0.5
62 RTN

[Dieter]

(12-18-2013 12:52 AM)Namir Wrote:  Bernoulli Numbers using Series Approximations

Hi Namir,

I think the RCL 00 in your line 31 should read RCL 01, right?

Then I tried my own implementation of this algorithm. I changed the way the tolerance is handled. If the user does not enter one, the current value in R01 is used. If this value is >=1 or <0, the tolerance defaults to 1E-5. Also the number of iterations is calculated beforehand so that the loop requires no tests and a simple DSE counter will do.

The following version uses only R00 - R02, and even R02 may be replaced with a stack register if the calculation of the factor starting at line 25 is placed behind the loop. Also no flags are required.

Code:


01 LBL"BN"
02 RCL 01
03 "TOL=?"
04 PROMPT
05 FRC
06 LASTX
07 X=Y?   ' tolerance >=1
08 X<=0?  ' or =< 0 ?
09 1E-5   ' set default tolerance
10 STO 01
11 "N=?"
12 PROMPT
13 STO 00
14 1
15 RCL 00
16 X=Y?
17 -,5
18 X<0?
19 RTN
20 2
21 MOD
22 -
23 X=0?
24 RTN
25 RCL 00
26 FACT
27 PI
28 ST+ X
29 RCL 00
30 Y^X
31 /
32 ST+ X
33 RCL 01
34 RCL 00
35 1/X
36 CHS
37 Y^X
38 INT
39 STO 02   ' set number of iterations
40 CLX
41 ISG 02
42 LBL 01
43 RCL 02
44 RCL 00
45 CHS
46 Y^X
47 +
48 DSE 02
49 GTO 01
50 *
51 RCL 00
52 -4
53 MOD
54 SIGN
55 *
56 CHS
57 END

What do you think? Any errors or things to change?

Dieter

[Dieter]

(12-18-2013 12:52 AM)Namir Wrote:  Bernoulli Numbers using Series Approximations

Hi again,

I could not resist and tried a new and improved version with a different approach. Since the summation loop will take substantial time for small n, all cases up to n = 8 are handled explicitely. For n >= 10 the result is evaluated iteratively. The accumulated terms add up to a value slightly higher than 1, so the smalles relevant term is ULP(1) = 1E-9 and the number of required terms is known. Even in the worst case (n = 10) not more than eight loops are required (8^-10 = 9,3E-10) and the result shows up after approx. 8 seconds.

If the machine value for Pi is sufficiently exact, the error should not exceed a few ULP. So there is no need to specify a tolerance. For the 10-digit 41-series where Pi is relatively inexact, a final correction is applied (see below).

Here is the algorithm in pseudocode:

Code:


if n =< 1 then return 1 - 1,5*n  ' n = 0 or 1: b = 1 resp. -1/2

if n is odd then return 0        ' n > 1 and odd: b = 0

if n =< 8 then                   ' n = 2, 4, 6, 8
   b = 1/ABS(3*(n-6)^2-42)       ' |b| = 1/6, 1/30, 1/42, 1/30
   goto setsign
                                 ' n>=10: use iteration
imax = CEIL(ULP(1)^(-1/n))       ' largest i where i^-n is
                                 ' significant in a sum near 1

for i = imax downto 1            ' improve accuracy by summing up
    sum = sum + i^(-n)           ' the smallest terms first
next i

b = sum * 2 * n! / (2*pi)^n

:setsign:
sign = sign41(n mod -4)   ' sign41 = 1 if argument is zero

return b * sign

And here is my current HP41 implementation:

Code:


01 LBL"BN2"
02 STO 00
03 1
04 RCL 00
05 X>Y?
06 GTO 00
07 1,5
08 *
09 -
10 GTO 04  ' exit
11 LBL 00
12 2
13 MOD
14 -       ' 1 - n mod 2
15 X=0?    ' = 0 if n is odd
16 GTO 04  ' return zero and exit
17 8
18 RCL 00
19 X>Y?
20 GTO 01  ' use iteration for n > 8 (i.e. n>=10)
21 6
22 -
23 X^2     ' use quadratic formula for n = 2, 4, 6, 8
24 3
25 *
26 42
27 -
28 ABS
29 1/X
30 GTO 03  ' adjust sign and exit
31 LBL 01
32 1E-9    ' use 1E-11 on 12-digit machines
33 RCL 00
34 CHS
35 1/X
36 Y^X
37 INT     ' largest i where i^(-n) >= 1E-9
38 0
39 ISG Y   ' start summation loop with next larger i
40 LBL 02
41 RCL Y
42 RCL 00
43 CHS
44 Y^X
45 +
46 DSE Y
47 GTO 02
48 RCL 00
49 FACT
50 PI
51 ST+ X   '  ten-digit value of pi has relative error 1,3 E-10
52 RCL 00
53 Y^X
54 /
55 ST+ X
56 *
57 RCL 00
58 1,3E-10  ' correct error in pi
59 *
60 RCL Y    ' on 12-digit machines
61 *        ' omit steps 57...62
62 +
63 LBL 03   ' adjust sign
64 RCL 00
65 -4
66 MOD
67 SIGN     ' returns 1 if n mod 4 = 0, and -1 otherwise
68 *        ' this is HP41-specific and will not work elsewhere
69 CHS
70 LBL 04   ' common exit point for all results
71 END

There is a slight tweak in line 57 ff. Since Pi does not round very well to ten digits, the value of pi has a relative error of 1,3E-10 that would show up in the last two digits. This error is compensated so that the final result should be correct within a few ULP. At least I did not notice any larger errors. ;-)

The domain for even n is limited to values up to 68, as 70! would exceed the 41's working range. B(68) is returned accurately as -2,625771029 E+42.

Dieter

[Namir]

Thanks for the correction and for your versions of the two listings.