02-09-2014, 04:59 PM
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}\).
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
\( 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