Bernoulli Numbers
12-18-2013, 12:52 AM (This post was last modified: 01-30-2014 09:24 PM by Namir.)
Post: #1 Namir Senior Member Posts: 823 Joined: Dec 2013
Bernoulli Numbers
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
01-28-2014, 09:21 PM
Post: #2
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Bernoulli Numbers
(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
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,397 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
01-30-2014, 09:34 PM
Post: #4 Namir Senior Member Posts: 823 Joined: Dec 2013
RE: Bernoulli Numbers
Thanks for the correction and for your versions of the two listings.
 « Next Oldest | Next Newest »