Bernoulli Numbers
01-30-2014, 02:37 PM (This post was last modified: 01-30-2014 02:56 PM by Dieter.)
Post: #3
 Dieter Senior Member Posts: 2,398 Joined: Dec 2013
RE: Bernoulli Numbers
(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
 « Next Oldest | Next Newest »

 Messages In This Thread Bernoulli Numbers - Namir - 12-18-2013, 12:52 AM RE: Bernoulli Numbers - Dieter - 01-28-2014, 09:21 PM RE: Bernoulli Numbers - Dieter - 01-30-2014 02:37 PM RE: Bernoulli Numbers - Namir - 01-30-2014, 09:34 PM

User(s) browsing this thread: 1 Guest(s)