Bernoulli numbers and large factorials
02-09-2014, 04:59 PM
Post: #1
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
Bernoulli numbers and large factorials
In December Namir posted a program that evaluates the Bernoulli numbers on the HP-41. It can be found in the HP-41C Software Library. In January I suggested another version based on the same approach:

$$B_n = 2 n! {(2\pi)^{-n}} \sum\limits_{i=1}^{\infty}i^{-n}$$

For large n the sum rapidly approaches 1 so that at some point it may be omitted. On the other hand, for n as low as 10 only a few terms are required for the usual 10 or 12 digit precision. $$B_{0...8}$$, which would require a substantial number of terms, may be given directly. So this is the easy part.

Both Namir's and my solution suffer from a limitation due to overflow in the factorial function. For HP's 10-digit calculators with working range up to 9,999...E99 the limit is 69!, while most 12-digit devices with their upper limit at 9.999...E499 will work up to 253!. However, this is significantly less than the largest possible n that returns a Bernoulli number within the working range. Here, results up to $$B_{116}$$ resp. $$B_{372}$$ would be possible.

There is an elegant way to overcome this problem if the calculator offers a permutation function (nPr). The essential idea is to split n! into two factors a and b which both fall within the calculator's working range. If factorials up to 253! are possible (the largest value below 9,999...E499) this can be done as follows:

\begin{align*} n! &= \frac{n!}{253!} · 253!\\ &= \frac{n!}{(n - (n - 253))!} · 253!\\ &= Perm(n, n - 253) · 253! \end{align*}

If n is less than 253, this constant is simply replaced with 0:

$$Perm(n, n - 0) · 0! = n! · 1 = n!$$

This method may also be written as follows:

let r = max(0, n - 253)
let a = Perm(n, r)
let b = (n-r)!

Then n! = a · b

This way the factor $$2 n! {(2\pi)^{-n}}$$ can be evaluated as $$2 a {(2\pi)^{-n} b}$$  to avoid overflow. On calculators with a limit at 9,999...E99, for instance the 15C, simply use 69 instead of 253.

Here is a complete program for the HP-35s. It evaluates all Bernoulli numbers within the working rage, i.e. from $$B_0$$ to $$B_{372}$$.

Code:
B001 LBL B B002 ABS B003 IP B004 STO N B005 1 B006 x>y? B007 RTN          ' n=0:  B = 1 B008 RCL N B009 x≠y? B010 GTO B015 B011 + B012 1/x B013 +/- B014 RTN          ' n=1:  B = -1/2 B015 2 B016 RMDR B017 - B018 x=0? B019 RTN          ' n is odd:  B = 0 B020 9 B021 RCL N B022 x>y? B023 GTO B035 B024 RCL N B025 6 B026 - B027 x² B028 3 B029 x B030 42 B031 - B032 ABS          '  n  = 2, 4, 6, 8: B033 1/x          ' |B| = 1/6, 1/30, 1/42, 1/6 B034 GTO B076 B035 RCL N        ' n ≥ 10:  use Bernoulli formula B036 253 B037 - B038 x<0? B039 CLx B040 nPr B041 RCL N B042 LASTx B043 - B044 ! B045 π B046 ENTER B047 + B048 RCL N B049 +/- B050 y^x B051 5E-14       ' adjust (2pi)^-n B052 RCLx N B053 x<>y B054 x B055 LASTx B056 + B057 x B058 x B059 ENTER B060 + B061 1E-12 B062 RCL N B063 +/- B064 x√y B065 IP B066 STO I       ' largest i where  i^(-n)  >  0,1 ULP B067 RCL- I      ' sum = 0 B068 RCL I B069 RCL N B070 +/- B071 y^x B072 + B073 DSE I B074 GTO B068 B075 x B076 RCL N    ' adjust sign B077 RCL N B078 4 B079 RMDR B080 x=0? B081 RCL- N B082 SGN B083 R↑ B084 x B085 RTN

The adjustment in line 051...056 tries to reduce the error in $$(2\pi)^{-n}$$. On the 35s I could not find errors larger than some units in the last place. Other calculators may require a different correction, or it may even be omitted completely if a somewhat larger error is acceptable.

Usage:
Enter n  [XEQ] B [ENTER]
=> display shows n and Bn.

Execution time:
within approx. 3 seconds (at n = 10).

Examples:

4 [XEQ] B [ENTER]
=> -0,333333333333

32 [XEQ] B [ENTER]
=> -15.116.315.767,1

100 [XEQ] B [ENTER]
=> 2,83322495708 E+78

372 [XEQ] B [ENTER]
=> -5,58475372908 E+499

Up to n = 26 the result may also be viewed in fraction mode. If no limitations are set (Flag 8 and 9 clear, /c ≥ 2730) the display shows the exact representation of the respective Bernoulli number:

22 [XEQ] B [ENTER]
=> 6.192,12318839

[FDISP] => 6192 17 / 138  ' exact result is  $$6192 \frac{17}{138}$$ or $$\frac{854513}{138}$$
[RND]       ' round to exact result
[FDISP] => 6.192,12318841

Dieter
02-09-2014, 07:47 PM
Post: #2
 Tugdual Senior Member Posts: 764 Joined: Dec 2013
RE: Bernoulli numbers and large factorials
Nice job
I checked the values on Wolfram Alpha and you are prettyu close (for the 10 first digits over 500 ;-) ).

The Prime has Bernouilli numbers calculated internally but unfortunately the Bn function is not exposed in the library. So we can use the formula Bernoulli(x) := –x*Zeta(1–x) (Thanks Joe Horn)

I tried it and couldn't calculate B372; the bigger number is B370.
02-09-2014, 08:01 PM
Post: #3
 Marcus von Cube Senior Member Posts: 760 Joined: Dec 2013
RE: Bernoulli numbers and large factorials
The WP 34S has Bn built-in. In double precision mode I can get B2122. Larger arguments return 0 (which can be considered a bug).

Marcus von Cube
Wehrheim, Germany
http://www.mvcsys.de
http://wp34s.sf.net
http://mvcsys.de/doc/basic-compare.html
02-09-2014, 08:42 PM
Post: #4
 Bunuel66 On Vacation Posts: 29 Joined: Jan 2014
RE: Bernoulli numbers and large factorials
Another approach would be to compute iteratively the factorial using 1/2π factor at every step. As the factorial as the same number of factor as the term (2π)^−n this would limit the overflow.

My two cents...
02-09-2014, 08:43 PM
Post: #5
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Bernoulli numbers and large factorials
(02-09-2014 07:47 PM)Tugdual Wrote:  Nice job
I checked the values on Wolfram Alpha and you are prettyu close (for the 10 first digits over 500 ;-) ).

Thank you. :-) I just compared all possible 188 non-zero results with the correctly rounded 12-digit values. If I got it right, more than 80% are within ±1 ULP and more than 95% within ±2 ULP. The rest (<5%) is off by ±3 or 4 ULP.

Dieter
02-09-2014, 09:26 PM (This post was last modified: 02-09-2014 10:13 PM by Dieter.)
Post: #6
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Bernoulli numbers and large factorials
(02-09-2014 08:01 PM)Marcus von Cube Wrote:  The WP 34S has Bn built-in. In double precision mode I can get B2122. Larger arguments return 0 (which can be considered a bug).

On my 34s (v. 3.2 3405) both B2124 and Zeta (-2123) still return a result, while beyond that a "+∞ Error" is displayed. Do you really get a zero here?

If 11 valid digits are sufficient, the 35s program does a good job and, compared to the 34s, it is really fast. I wonder how the 34s will perform with the same algorithm in user code. OK, 34 digits for n as low as 10 or 12 will take somewhat longer. ;-)

EDIT: The reason for the 34s limit at B2122 probably is the same as the one mentioned in my original post: it's the factorial function. In DP mode the 34s still can handle 2122! but 2124! will cause an overflow. This could be overcome by using the permutation function. A quick-and-dirty test confirmed that B2776 = 1,01268...E+6140 can be done. The result is returned in about a second.

Dieter
02-09-2014, 09:44 PM
Post: #7
 Paul Dale Senior Member Posts: 1,770 Joined: Dec 2013
RE: Bernoulli numbers and large factorials
The 34S is computing the Bernoulli numbers from the zeta function. A short series expansion like you've got here will be faster I expect.

- Pauli
02-12-2014, 08:56 PM
Post: #8
 Dieter Senior Member Posts: 2,397 Joined: Dec 2013
RE: Bernoulli numbers and large factorials
(02-09-2014 09:44 PM)Paul Dale Wrote:  The 34S is computing the Bernoulli numbers from the zeta function. A short series expansion like you've got here will be faster I expect.

Well, this series expansion, i.e. the sum in the formula, actually is the Zeta function. ;-)

\begin{align*} B_n &= 2 n! {(2\pi)^{-n}} \sum\limits_{i=1}^{\infty}i^{-n}\\ &= 2 n! {(2\pi)^{-n}} \zeta (n) \end{align*}

The 35s program works so fast because the larger n gets, the less terms are required. The following table shows the number of terms needed for an error of at most 0,1 ULP in Zeta:

Code:
 n      10      12      16    34 digits ---------------------------------------  8      17      31     100     17782 10      10      15      39      2511 12       6      10      21       681 14       5       7      13       268 16       4       5      10       133 20       3       3       6        50 24       2       3       4        26 30       2       2       3        13 40       1       1       2         7

This also explains why up to n = 8 the result is given directly. Otherwise the number of required terms would increase rapidly.

I tried a program with the same algorithm on the 34s. In SP mode it is much faster than the internal Bernoulli function. For large n the result appears within a fraction of a second. As you will expect after a look at the above table, DP mode with 34 digit precision is a different story, at least for small n. ;-)

Dieter
 « Next Oldest | Next Newest »