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
 « Next Oldest | Next Newest »

 Messages In This Thread Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-09-2014 07:12 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-10-2014, 12:20 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-10-2014, 08:23 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-12-2014, 04:34 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-14-2014, 10:05 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-14-2014, 09:37 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-14-2014, 10:00 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-14-2014, 10:21 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-14-2014, 10:04 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-15-2014, 09:32 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-15-2014, 12:26 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-15-2014, 06:12 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Thomas Klemm - 03-15-2014, 07:17 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-15-2014, 12:56 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-15-2014, 08:47 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-18-2014, 09:24 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-19-2014, 06:47 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-19-2014, 01:37 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-19-2014, 07:32 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-19-2014, 08:31 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-20-2014, 06:46 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-20-2014, 12:44 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-20-2014, 02:21 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-20-2014, 09:07 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Ángel Martin - 03-21-2014, 05:54 AM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-21-2014, 01:37 PM RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"? - Dieter - 03-28-2014, 07:56 PM

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