Bernoulli numbers and large factorials

02092014, 04:59 PM
Post: #1




Bernoulli numbers and large factorials
In December Namir posted a program that evaluates the Bernoulli numbers on the HP41. It can be found in the HP41C 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 10digit calculators with working range up to 9,999...E99 the limit is 69!, while most 12digit 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 = (nr)! 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 HP35s. It evaluates all Bernoulli numbers within the working rage, i.e. from \(B_0\) to \(B_{372}\). Code: B001 LBL B 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 

02092014, 07:47 PM
Post: #2




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. 

02092014, 08:01 PM
Post: #3




RE: Bernoulli numbers and large factorials
The WP 34S has Bn builtin. 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/basiccompare.html 

02092014, 08:42 PM
Post: #4




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... 

02092014, 08:43 PM
Post: #5




RE: Bernoulli numbers and large factorials
(02092014 07:47 PM)Tugdual Wrote: Nice job Thank you. :) I just compared all possible 188 nonzero results with the correctly rounded 12digit 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 

02092014, 09:26 PM
(This post was last modified: 02092014 10:13 PM by Dieter.)
Post: #6




RE: Bernoulli numbers and large factorials
(02092014 08:01 PM)Marcus von Cube Wrote: The WP 34S has Bn builtin. 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 quickanddirty test confirmed that B2776 = 1,01268...E+6140 can be done. The result is returned in about a second. Dieter 

02092014, 09:44 PM
Post: #7




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 

02122014, 08:56 PM
Post: #8




RE: Bernoulli numbers and large factorials
(02092014 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 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 »

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