Accurate Bernoulli numbers on the 41C, or "how close can you get"?
03-09-2014, 07:12 PM
Post: #1
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
Accurate Bernoulli numbers on the 41C, or "how close can you get"?
In the last weeks there have been some discussions regarding various ways of determining Bernoulli numbers on the 41-series and other calculators. The usual formulas included powers with exponents greater than 100, leading to reduced accuracy since an exact result would require at least twelve or thirteen digits for the base as opposed to the ten we have. Another problem is the available working range, so that the used algorithm has to make sure no intermediate result exceeds the limit at 9,999...E99.

So I wondered if there might be a way of evaluating all possible Bernoulli numbers within the working range sufficiently fast and, more important, as accurately as possible. Which leads to the question: how close can you get in the face of accumulating roundoff errors? Even a simple multiplication can be surprisingly inaccurate, the result may be off by up to 5 units in the last place. Try a simple $$\pi·\pi$$ or $$e·e$$, and the result on a correctly working 10-digit calculator is 3 ULP high or low. So far, so bad.

Here is the approach used in the following program. As usual, the lower Bernoulli numbers B0 to B8 are given directly. For n = 2...8 a simple quadratic equation can do the trick. For n = 10 to 116 (largest value within the 41's working range) the following formula was used:

$$\large B_n = 4 \pi · \zeta(n) · e^{(0,5 + \frac{1}{12n} - \frac{1}{360n^3} + \frac{1}{1260n^5} + ...)} · (\frac{n}{2 \pi e})^{n+0,5}$$

For n ≥ 10 and 10-digit accuracy three terms of the series in the exponent of the e-function are sufficient.
The expression $$\frac{1}{12n} - \frac{1}{360n^3} + \frac{1}{1260n^5}$$ can be evaluated as $$\frac{210n^4 - 7n^2 + 2}{2520n^5}$$.

A literal implementation of the complete formula would yield results with substantial errors. At least the last two digits would be off. So a different way to handle this formula had to be found.

Within the relevant domain, the factors at the left ($$4\pi=10·0,4\pi, \zeta$$ and the exponential function) all start with 1. They do not vary much for n = 10...116:

$$B_n = 10 · 1,256... · 1,000... · 1,65... · (\frac{n}{2 \pi e})^{n+0,5}$$

The basic idea now is to evaluate all three factors minus one so that one additional digit is gained. Obtaining $$\zeta - 1$$ is trivial, and for the exponential function there is a dedicated $$e^x-1$$ command. The multiplication of three values close to 1 can be done in a way that preserves one additional digit of working precision. Since the product of the three factors is something between 2,07 and 2,09, the program even tries to calculate half of this minus 1 (and finally multiplies this +1 with twice the power), so that again a precious digit is saved. The program uses a 9-digit approximation of $$0,4\pi - 1 = rad(72°)-1$$ which is slightly low, so a correction term is applied. Its exact value should be near 7,2E-10, but tests showed that in this case even better accuracy is obtained with a slighty lower value close to 6E-10 (cf. line 89).

Now let's look at the power at the right. For a correct 10-digit result, the base would have to carry at least 12 or 13 digits. Here is how this is accomplished in the program:

$$(\frac{n}{2 \pi e})^{n+0,5}$$
$$= (n · 0,05854983152432)^{n+0,5}$$
$$\approx (n · 0,05854983)^{n+0,5} + 1,52432·10^{-9}·n·(n+0,5)·(n · 0,05854983)^{n-0,5}$$

For n = 10...116 the base of the first power carries at most 9 digits, so both the base and the exponent are exact. However, the 41's power function sometimes truncates its result instead of rounding it, so the constant 1,52432E-9 is better rounded up to 1,5244E-9.

Here is the 41C code:
Code:
 01  LBL"BN"  02  ABS  03  INT  04  STO 00  05  SIGN  06  RCL 00  07  X>Y?  08  GTO 01  09  1,5  10  *  11  -  12  GTO 99  13  LBL 01  14  2  15  MOD  16  -  17  X=0?  18  GTO 99  19  9  20  RCL 00  21  X>Y?  22  GTO 02  23  6  24  -  25  X^2  26  3  27  *  28  42  29  -  30  ABS  31  1/X  32  GTO 98  33  LBL 02  34  5 E-10  35  RCL 00  36  CHS  37  1/X  38  Y^X  39  INT  40  STO 01  41  0  42  ISG Y  43  LBL 03  44  RCL Y  45  RCL 00  46  CHS  47  Y^X  48  +  49  DSE Y  50  DSE 01  51  GTO 03  52  RCL 00  53  X^2  54  ENTER  55  ENTER  56  210  57  *  58  7  59  -  60  *  61  2  62  +  63  2520  64  /  65  RCL 00  66  ENTER  67  X^2  68  X^2  69  *  70  /  71  ,5  72  +  73  E^X-1  74  ST* Z  75  +  76  +  77  ENTER  78  ENTER  79  1  80  -  81  72  82  D-R  83  FRC  84  +  85  X<>Y  86  LASTX  87  *  88  +  89  5,8 E-10  90  +  91  2  92  /  93  STO 01  94  RCL 00  95  ,05854983  96  *  97  ENTER  98  ENTER  99  RCL 00 100  ,5 101  + 102  Y^X 103  ST+ X 104  ENTER 105  ENTER 106  R^ 107  / 108  1,5244 E-9 109  * 110  RCL 00 111  2 112  / 113  RCL 00 114  X^2 115  + 116  * 117  + 118  RCL 01 119  X<>Y 120  * 121  LASTX 122  + 123  10 124  * 125  LBL 98 126  RCL 00 127  -4 128  MOD 129  SIGN 130  * 131  CHS 132  LBL 99 133  END

One may now ask if the result is worth all the effort. I think it is. In total there are 60 possible non-zero results within the 41's working range (n = 0, 1, 2, 4, 6, 8, ..., 114, 116). The program returns 45 of these correctly rounded or truncated after 10 digits. The rest is 1 ULP high or low. I did not find any larger errors. In other words: the results are close to machine accuracy.

BTW, while the largest possible result is B116, the program can also provide B118. The expected OUT OF RANGE error appears in the last calculation step when the program tries to multiply X by 10. At this point, pressing [X<>Y] reveals B118 as 6,116052000E+100. ;-)

Of course suggestions for improvements are always welcome.

Dieter
03-10-2014, 12:20 AM
Post: #2
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
Did you ever consider to write the program in MCODE? That provides the possibility to use 13 digits for your calculations and only round the result to 10 digits.

Cheers
Thomas
03-10-2014, 08:23 PM
Post: #3
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-10-2014 12:20 AM)Thomas Klemm Wrote:  Did you ever consider to write the program in MCODE? That provides the possibility to use 13 digits for your calculations and only round the result to 10 digits.

Getting 10 valid digits out of 13 digits working precision is not quite a challenge. ;-)

The basic question here was: how close can you get without any additional guard digits. Otherwise things are trivial: I tried the basic formula on a 12-digit 35s. Even with sloppy programming 10 valid digits ±1 ULP are not very hard to achieve. So with 13 digit precision things become trivial if the goal is merely a correct 10-digit result.

But nevertheless I would love to try MCODE programming. Years ago I always enjoyed combining efficient x86 assembler code with the beauty of Pascal. But honestly, I do not know the slightest thing about MCODE. I do not even how what is required to build a working programming environment. Maybe there is some basic info somewhere to start with?!

BTW, in the meantime there is a shorter and even slightly faster Bernoulli program with comparable accuracy: ±1 ULP over the whole working range. Yes, with a standard 10-digit 41C. ;-)

Dieter
03-12-2014, 04:34 PM
Post: #4
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-10-2014 08:23 PM)Dieter Wrote:  Getting 10 valid digits out of 13 digits working precision is not quite a challenge. ;-)

(03-12-2014 01:45 AM)Paul Dale Wrote:  Having extra guard digits during a computation is a godsend.

But nevertheless I would love to try MCODE programming.

You'll find all you need on Warren Furlow's HP-41 web-site.

Have fun
Thomas
03-14-2014, 10:05 AM
Post: #5
 Ángel Martin Senior Member Posts: 1,261 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
It's slightly more involved than just applying a cookie-cutter recipy, even if it's all published somewhere (like TOS). The 13-digit routines allow for intermediate calculations done with more precision, so the final result is more accurate (or less inaccurate, in a more strict speak). There are many 13-digit versions of the standard math routines in the 41 OS, a good starting point could be the mini-paper I prepared a while ago; available at TOS (can't paste the link here as you probably know).

I've attached it to this post as well (hope it works)

Cheers,
ÁM

Attached File(s)
03-14-2014, 09:37 PM
Post: #6
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-12-2014 04:34 PM)Thomas Klemm Wrote:  You'll find all you need on Warren Furlow's HP-41 web-site.
(03-14-2014 10:05 AM)Ángel Martin Wrote:  It's slightly more involved than just applying a cookie-cutter recipy, even if it's all published somewhere (like TOS).

Yes, I am absolutely sure it is all "somewhere". ;-) But what I am looking for is something like a structured introduction to MCODE programming: what are the hardware options, what else is required, what is covered by the instruction set and what are basic programming techniques. I am quite familiar with FOCAL programming, and I also appreciate the the benefits of doing things in machine code, so MCODE sounds like an interesting new project to me. However, I need some basic info first.

Quote:The 13-digit routines allow for intermediate calculations done with more precision, so the final result is more accurate (or less inaccurate, in a more strict speak). There are many 13-digit versions of the standard math routines in the 41 OS, a good starting point could be the mini-paper I prepared a while ago; available at TOS (can't paste the link here as you probably know).
I've attached it to this post as well (hope it works)

Thank you very much, and yes, it worked. That's exactly what I am looking for. In a way I am a kind of "accuracy junkie" - that's why I love the 34s project and the option of 34-digit precision for certain things. Defining special functions like the Normal distribution quantile, Lambert's W or just the Bernoulli numbers will full 10-digit accuracy is virtually impossible in FOCAL. With 13 digits all this becomes easy.

Now - where and how to start?

Dieter

By the way, do you know of a kind of overview with the actual accuracy of the 41's internal 13-digit routines?
03-14-2014, 10:00 PM
Post: #7
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-14-2014 09:37 PM)Dieter Wrote:  Now - where and how to start?

Software Development Kit
For the HP-41 Release 6
PROGRAMMER'S MANUAL

cf. APPENDIX I – Step by Step Example

HTH
Thomas
03-14-2014, 10:04 PM
Post: #8
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-14-2014 10:05 AM)Ángel Martin Wrote:  available at TOS (can't paste the link here as you probably know).
Not a problem: HP41 OS 13-digit Math Routines

Thanks again for the compilation!

Cheers
Thomas
03-14-2014, 10:21 PM
Post: #9
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-14-2014 09:37 PM)Dieter Wrote:  By the way, do you know of a kind of overview with the actual accuracy of the 41's internal 13-digit routines?

Probably not what you want:
HEWLETT-PACKARD JOURNAL
cf. Algorithms and Accuracy in the HP-35

But I just remembered that picture of the saw-tooth diagram.

Cheers
Thomas
03-15-2014, 09:32 AM
Post: #10
 Ángel Martin Senior Member Posts: 1,261 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
There's an interesting article (written by Dennis W. Harms) that broadly describes the accuracy of the math algorithms on the HP-67 - which in my understanding were pretty much the same ones also used by the 41C. You can read it in the November 1976 issue of the HP Journal, pages 16 and 17.

Cheers,
'AM
03-15-2014, 12:26 PM
Post: #11
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-15-2014 09:32 AM)Ángel Martin Wrote:  You can read it in the November 1976 issue of the HP Journal, pages 16 and 17.
Original article in Advanced Pocket Calculators

Dennis W. Harms, pictured 30 years later:
03-15-2014, 12:56 PM
Post: #12
 Ángel Martin Senior Member Posts: 1,261 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
Thomas, thanks for the links, I didn't have the original article.
Your historian archivist skills largely exceed mine
03-15-2014, 06:12 PM
Post: #13
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-15-2014 12:26 PM)Thomas Klemm Wrote:
(03-15-2014 09:32 AM)Ángel Martin Wrote:  You can read it in the November 1976 issue of the HP Journal, pages 16 and 17.
Original article in Advanced Pocket Calculators

Thank you very much, especially for the "Advanced Pocket Calculators" version. I read the HP journal version of this article many years ago, but it does not answer the question regarding the accuracy of the 13-digit internal functions of the HP-41, -67 or any other calculator. What the article essentially says is that 13 digits (HP-91 and newer) allow (nearly) correct 10-digit results, which had errors before when no guard digits were used. Or simply "13 digits give better results than 10 digits". Which nobody will deny. ;-)

Dieter
03-15-2014, 07:17 PM
Post: #14
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-15-2014 06:12 PM)Dieter Wrote:  but it does not answer the question regarding the accuracy of the 13-digit internal functions of the HP-41
There's most probably no information available concerning this question. Thus there is much to do. I suggest to use the trace which is available in Nonpareil to figure out the 13-digit internal results of the operations you are interested in. You probably have some corner cases already in mind such as cos(x) close to $$\frac{\pi}{2}$$ or ln(x) close to 1. Or maybe the SDK41 provides a decent debugger? Let me know if you need detailed instructions. IIRC Nonpareil doesn't support tracing out of the box. You have to compile it setting the option has_debugger_gui.

Cheers
Thomas
03-15-2014, 08:47 PM (This post was last modified: 03-15-2014 08:50 PM by Ángel Martin.)
Post: #15
 Ángel Martin Senior Member Posts: 1,261 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
V41 does have a tracer and breakpoint capability. I do lots of troubleshooting on it, and it can also be used to see the intermediate results - provided you have the source code of course. So it's a bit of a catch-22, but once you're in the game it's very rewarding.
03-18-2014, 09:24 PM
Post: #16
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
Thomas and Ángel: Thank you very much for some first information on how to start with MCODE. It seems to be a long way to go. On the other hand I just read through the Sandmath manual, and all this looks very promising.

BTW, Ángel: Sandmath seems to lack a Normal Quantile function. I do not think it's as complicated as Appendix 9 suggests: on the 34s it can be done with 30+ digit accuracy in merely two iterations. At the moment I am playing around a bit with V41. If there is a way to have Sandmath's ERF in V41 I could see what can be done in user code. ;-) How fast does ERF execute on a real 41?

Dieter
03-19-2014, 06:47 AM (This post was last modified: 03-19-2014 06:52 AM by Ángel Martin.)
Post: #17
 Ángel Martin Senior Member Posts: 1,261 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
Hi Dieter, glad to see you find the SandMath worth looking at. Which revision are you using? Are there two appendices 9a and 9b, or only one appendix 9?

http://hp41.claughan.com/file/SANDMATH_4...Manual.pdf

I say this because in the latest versions (Revision "M", see above link) there are a couple of functions directly related to the quantile, one (ICPF) using the inverse error function on a general Normal distribution (s,m), and another (QNTL) using Halley's iterative method on a standard Normal (0,1). The functions are a bit buried in the -FACTORIAL group, in the sub-functions FAT

The manual has room for improvement, but isn't the quantile function QNTL what you're asking about? It is briefly (and poorly) documented in page 48. Below you can see the FOCAL code for this function, which uses Halley's method to solve for a CPF (Cumulative Probability Function) equation equal to the quantile's value. The convergence criteria is a poor E-7, probably not what you're looking for I'm afraid.

Code:
 01    LBL "QNTL" 02    STO 03 03    CLX 04    STO 04 05    LBL 01 06    0 07    ENTER^ 08    1 09    RCL 04 10    CPF 11    RCL 03 12    - 13    STO 05 14    0 15    ENTER^ 16    1 17    RCL 04 18    PDF 19    RCL 05 20    X<>Y 21    ST* Y 22    X^2 23    LASTX 24    RCL 04 25    * 26    RCL 05 27    * 28    2 29     / 30    + 31     / 32    ST- 04 33    ABS 34    E-7 35    X<Y? 36    GTO 01 37    END

As per ICPF, it's a much simpler code since IERF does all the work in there:

Code:
 01    LBL "ICPF"  02    ST+ X 03    E 04    - 05    IERF 06    2 07    SQRT 08    * 09    * 10    + 11    END

I use the CUDA library algorithms to calculate the inverse error function - it's a very fast implementation (yes, even on a normal HP-41) that takes a huge MCODE listing code.

I may be confusing subjects, statistics is not my forte (I'm really a jack of all trades and master of none, as I'm sure you have noticed ;-) I think this function is equivalent to the 38E's, but haven't looked at the 34S at all (don't have one myself).

l hope this clarifies, but let me know if you see any issues or would like additions/changes made.

Best,
'AM
03-19-2014, 01:37 PM (This post was last modified: 03-19-2014 06:54 PM by Dieter.)
Post: #18
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-19-2014 06:47 AM)Ángel Martin Wrote:  Hi Dieter, glad to see you find the SandMath worth looking at. Which revision are you using? Are there two appendices 9a and 9b, or only one appendix 9?

The version I had was 4.44.44, now updated to 5.55.55 that indeed has additional information on the Normal quantile function.

Actually I do not use Sandmath. I heard about it quite often here, and when I tried a first look at MCODE I came across the manual. At the moment this is all what I got. ;-)

(03-19-2014 06:47 AM)Ángel Martin Wrote:  I say this because in the latest versions (Revision "M", see above link) there are a couple of functions directly related to the quantile, one (ICPF) using the inverse error function on a general Normal distribution (s,m), and another (QNTL) using Halley's iterative method on a standard Normal (0,1).

Halley's method is nice, but it can be done better and faster. ;-)

I now remember that some time ago I set up an algorithm specially taylored for the HP41 and other 10-digit machines. It first calculates an initial estimate with a very simple (1, 2) rational function, requiring just four constants with few digits, and then one single correction step is sufficient for a result with error less than 0,055 ULP, i.e. 5,5 units in the 12th significant digit. All it needs is a sufficiently accurate CDF (cumulative distribution function).

The following method assumes p < 0,5. The rest is handled by CDF(-x) = 1 - CDF(x).

$$u := \sqrt{-2 ln p}$$

$$x := max( 0 , u - \large \frac{2,374 + 0,39 u}{1 + 1,11 u + 0,07157 u^2} )$$

$$t := \large \frac{CDF(x) - p}{PDF(x)}$$

Now apply a third order (!) correction, cf. Abramowitz & Stegun, 10th ed., p. 954:

$$x := x + t + \frac{x}{2}t^2 + \frac{2x^2+1}{6}t^3$$

If evaluated exactly, the result has an error within approx. +0 and -5,5 units in the 12th significant digit throughout the 41's working range (even down to 1E-121).

Here is an implementation in Visual Basic for Excel:
Code:
Function icdf10(p)   signflag = (p > 0.5)   If signflag Then q = 1 - p Else q = p   u = Sqr(-2 * Log(q))   x = u - (2.374 + 0.39 * u) / ((0.07157 * u + 1.11) * u + 1)   If x < 0 Then x = 0   t = (normcdf(x) - q) / normpdf(x)   x = t * t * t * (2 * x * x + 1) / 6 + t * t * x / 2 + t + x   If signflag Then x = -x   icdf10 = x End Function
Here normcdf() and normpdf() refer to the Normal CDF resp. PDF. Please note that the PDF refers to the upper integral, i.e. from x to infinity. Otherwise two small adjustments are required.

The essential magic is in the combination of a dedicated initial guess and a very effective correction afterwards. If the term at t³ is omitted, the result quite exactly matches that of Halley's method and the error is less than 2 units in the 9th digit.

All this should work for p = 0,2 as well as p = 1E-12 or 1E-99 (that's why $$erf^{-1}(2p-1)$$ is not a good idea). Cases close to the distribution center and and the lower working limit (1E-99) require some care at the calculation of CDF(x) - p since this difference may become numerically zero due to underflow or digit cancellation. But there are ways to handle this.

Dieter

Edit: corrected an error and tweaked one of the constants for even better accuracy.
03-19-2014, 07:32 PM
Post: #19
 Ángel Martin Senior Member Posts: 1,261 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
Will look at your method as soon as I get a chance, thanks for sharing it.

I only included the Halley's method QNTL for comparison reasons; the preferred approach is no doubt to use IPFC (based on the inverse error function) - much faster (really several orders of magnitude) and generally more accurate.

BTW you can use the SandMath on V41 when you want, the MOD file is available at TOS.

Cheers,
'AM
03-19-2014, 08:31 PM
Post: #20
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03-19-2014 07:32 PM)Ángel Martin Wrote:  ...the preferred approach is no doubt to use IPFC (based on the inverse error function) - much faster (really several orders of magnitude) and generally more accurate.

"Fast and accurate" always sounds very promising. ;-) Do you have any figures here? How fast is, for instance, qf(0,01) on an actual HP41? You say you used the inverse error function as implemented in the CUDA library. Do you have a link to this algorithm? I found some information on the library, but not on the actual algorithms.

Just for comparison: I currently use a method similar to the one posted above in combination with an optimized method for the CDF. In a regular FOCAL program, the CDF runs in about 3 - 12 seconds, the quantile function requires roughly 3 seconds more. Maybe you can estimate how fast this would run in MCODE.

Dieter
 « Next Oldest | Next Newest »

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