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,401 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)