(12C+) Bernoulli Number

[Gamo]

In need of the Bernoulli Number using HP-12C ?

Here is an attempt to generate a Bernoulli Number constant using 12C

Without a Pi function this program use 355/113 which give out about

4 to 5 digits precision.

-------------------------------------------------
To run:
If you need to know B10 divide it by 2 is 5

5 [R/S] display 0.07576 [R/S] 5 [X<>Y] 66

Answer: B10 is 0.07576 or in fraction is 5/66
--------------------------------------------------
B12

12 ÷ 2 = 6

6 [R/S] display 0.25311 [R/S] 61 [X<>Y] 241

Answer:
B12 since 12 is divisible by 4 answer is Negative

-0.25311 in fraction is -61/241
--------------------------------------------------
Remark:
To find B(n) divide it by 2 and calculate.
This program do not give answer of the alternate negative value
such as B2 = 1/6 where B4 = -1/30
For B(n) that divisible by 4 answer is "Negative"
--------------------------------------------------
Program:

Code:


2
x
STO 2
355
ENTER
113
÷
STO 3
1
STO 0
STO 1
-----------------
RCL 0  // line 16
2 x 1 +
RCL 2
CHS
Y^X
RCL 1
+
RCL 1
X<>Y
X≤Y
GTO 34
STO 1
1
STO+0
GTO 16
RCL 2  // Line 34
n!
2 x
RCL 1
x 2
RCL 2
Y^X
1 -
RCL 2
RCL 3
X<>Y
Y^X
x ÷  // B(n) constant end here
----------------
R/S  // Line 51
STO 0  // Decimal to Fraction start here
STO 1
0
STO 2
1
RCL 0
INTG  // Line 58
RCL 2
x +
STO 2
RCL 0
x
. 5 // decimal and five // Line 66
+
INTG
STO 3
RCL 2
÷
RND
RCL 0
RND
-  // Subtract sign
X=0
GTO 86
Rv  // Roll Down
RCL 2
X<>Y
RCL 1
FRAC
1/x
STO 1
GTO 58
RCL 2  // Line 86
RCL 3
GTO 00  // Line 88

Formula use to calculate Bernoulli Number

B(n) = [2(2n)! ÷ ((2^2n) - 1)(Pi^2n)] [1 + (1/3^2n) + (1/5^2n) + ...]

Gamo

[Albert Chan]

(07-27-2019 06:41 AM)Gamo Wrote:  Formula use to calculate Bernoulli Number

B(n) = [2(2n)! ÷ ((2^2n) - 1)(Pi^2n)] [1 + (1/3^2n) + (1/5^2n) + ...]

Hi, Gamo

I think there is a typo: B(n) should be abs(B(2n))

Can you provide source ?

From https://wstein.org/edu/fall05/168/projec...rnproj.pdf

zeta(x) = sum(k^-x, where k = 1 to ∞) = 1 / product(1 - p^-x, where p is prime)

B(2n) = (-1)^(n+1) * 2*(2n)!/(2 pi)^(2n) * zeta(2n)

[Gamo]

Hello Albert Chan

Thanks for the review.
Here is the formula:

[attachment=7559]

Look like this formula is only good for
Small B(n)

Gamo

[John Keith]

It looks like an exact formula to me but it is an infinite series, and I have no idea how fast it converges. The Wikipedia page on Bernoulli numbers has lots of useful information including some exact algorithms. The fastest methods do seem to require a lot of memory so they may not be suitable for a small calculator like the 12C.

[Albert Chan]

I lookup "A Source Book in Mathematics", chapter "On the Bernoulli numbers":

B(n) is just sum of k^n formula linear term coefficient.

Example, this is how B(6) is calculated, by doing k^6 forward difference
(see thread: https://www.hpmuseum.org/forum/thread-12...#pid110972)

1 64 729 4096 15625 46656 117649 // value of 1^6 to 7^6
63 665 3367 11529 31031 70993     // forward differences
602 2702 8162 19502 39962
2100 5460 11340 20460
3360 5880 9120
2520 3240
720

Sum of k^6 formula = \(1\binom{n}{1}+63\binom{n}{2}+602\binom{n}{3}+2100\binom{n}{4}+3360\binom{n}{5}​+2520\binom{n}{6}+720\binom{n}{7}\)

B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42

Correction:
B(n) is sum of k^n formula, Σ(i^k, i = 0 to n-1) linear term coefficient
B(1) = linear coefficient of n(n+1)/2 - n = 1/2 - 1 = -1/2
(see http://www.mikeraugh.org/Talks/Bernoulli...n-LACC.pdf, slide 31 to 34)

The Python code (next post), work for n=1 because it flipped sign for all odd n.
The "bug" actually fix the sign

[Albert Chan]

Python code for Bernoulli Number (also return the forward difference table)

Code:

from fractions import Fraction

def B_fdiff(n):
    lst = [1,1] + [i**n for i in range(2, n+2)]
    for i in range(n):
        for j in range(n,i,-1): lst[j+1] -= lst[j]  # forward diff
        lst[0] = Fraction(lst[i+2],i+2) - lst[0]    # B(n) value
    return lst

>>> B_fdiff(6)
[Fraction(1, 42), 1, 63, 602, 2100, 3360, 2520, 720]

>>> for i in range(13): print i, B_fdiff(i)[0]
...
0 1
1 -1/2
2 1/6
3 0
4 -1/30
5 0
6 1/42
7 0
8 -1/30
9 0
10 5/66
11 0
12 -691/2730

[Gamo]

I use this same program on RPN-67 app the HP-67 emulator.

This RPN-67 got the printer feature so I did the loop print out.

The result give good approximate answer from B2 to B112

B114 just give out incorrect answer to a whole number 107407655 the correct
answer should be 5.17567436175.....

Correct answer continue from B116 to B124

Then B126 and over all result are 1.000000000 ***
Is this have to do with the limited value of the Factorial function?

Clip: https://youtu.be/LH7yDHan08M

Gamo

[John Keith]

(07-28-2019 12:02 AM)Albert Chan Wrote:  I lookup "A Source Book in Mathematics", chapter "On the Bernoulli numbers":

B(n) is just sum of k^n formula linear term coefficient.

Example, this is how B(6) is calculated, by doing k^6 forward difference
(see thread: https://www.hpmuseum.org/forum/thread-12...#pid110972)

1 64 729 4096 15625 46656 117649 // value of 1^6 to 7^6
63 665 3367 11529 31031 70993     // forward differences
602 2702 8162 19502 39962
2100 5460 11340 20460
3360 5880 9120
2520 3240
720

Sum of k^6 formula = \(1\binom{n}{1}+63\binom{n}{2}+602\binom{n}{3}+2100\binom{n}{4}+3360\binom{n}{5}​+2520\binom{n}{6}+720\binom{n}{7}\)

B(6) = Linear term coefficient = 1/1 - 63/2 + 602/3 - 2100/4 + 3360/5 - 2520/6 + 720/7 = 1/42

That is a very neat method, I was not aware of that one. However, it seems that all similar exact methods require n+(n-1) storage registers to calculate B(n) since one needs to keep the (n-1)th row of the difference table in memory while calculating the nth row.

EDIT: I tried your program as well as the Akiyama-Tanagiwa method as used in the third program here on the HP-48G and both methods fail due to catastrophic rounding error for n>10. These methods may only be practical for languages that use exact rational arithmetic.

[Albert Chan]

Knowing the pattern of s_k(n) = (Σi^k, k=0 to n-1) = n^(k+1)/(k+1) - n^k/2 + k/12 * n^(k-1) + 0 * x^(k-2) + ...,
we can get Bernoulli constants another way.

Example: to get upto B(6), s₆(n) = n^7/7 - n^6/2 + n^5/2 + a*n^3 + b*n, where a, b are unknown

→ a*n^2 + b = s₆(n)/n - n^4*(n^2/7 - (n-1)/2)

For n=1, eqn1 = a + b = 0/1 - 1*(1/7 - 0/2) = -1/7
For n=2, eqn2 = 4a + b = 1/2 - 16*(4/7 - 1/2) = -9/14

eqn2 - eqn1: 3a = -9/14 + 2/14 = -1/2 → a = -1/6

B(6) = b = -1/7 - a = 1/42
B(4) = a * 3 / \(\binom{6}{4}\) = -1/30
B(2) = (1/2) * 5 / \(\binom{6}{2}\) = 1/6
B(1) = (-1/2) * 6 / \(\binom{6}{1}\) = -1/2
B(0) = (1/7) * 7 / \(\binom{6}{0}\) = 1