Bernoulli Numbers - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: HP-41C Software Library (/forum-11.html) +--- Thread: Bernoulli Numbers (/thread-143.html) |
Bernoulli Numbers - Namir - 12-18-2013 12:52 AM Bernoulli Numbers using Series Approximations Algorithm Given integer n and tolerance toler Code: IF n=1 THEN Return -0.5 HP-41C Implementation Memory Map R00 = n R01 = toler R02 = i R03 = Sum R04 = term Implementation Code: 1 LBL "BERSER" RE: Bernoulli Numbers - Dieter - 01-28-2014 09:21 PM (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:
What do you think? Any errors or things to change? Dieter RE: Bernoulli Numbers - Dieter - 01-30-2014 02:37 PM (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:
And here is my current HP41 implementation: Code:
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 RE: Bernoulli Numbers - Namir - 01-30-2014 09:34 PM Thanks for the correction and for your versions of the two listings. |