(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-Tanigawa 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

[Albert Chan]

(07-28-2019 11:21 AM)John Keith Wrote:  EDIT: I tried your program as well as the Akiyama-Tanigawa method as used in the third program here on the HP-48G and both methods fail due to catastrophic rounding error for n>10.

If we use scaled numbers, we can push it a bit more.
Example, this is Akiyama-Tanigawa method, scaled by denominator LCM.

Note: XCas 1.9.0 default float now is 48 bits, truncated rounding.
Note: round(L[0]) does not change L[0] value, just its type, real → integer.

XCas> m := 22;
XCas> h := lcm(range(2, m+2))         → 5354228880
XCas> L := float(h ./ range(2, m+2))
XCas> for(k:=2; k<=m; k++) {
                L := -deltalist(L) .* range(1,len(L));
                if (even(k)) print(k, round(L[0])/h);
          };

2,      1/6
4,     -1/30
6,      1/42
8,     -1/30
10,     5/66
12,    -691/2730
14,     7/6
16,    -3617/510
18,     43867/798
20,    -174611/330
22,     854513/138

48 bits float mangaged to get upto B(22)
IEEE double can reach B(24), all exact fractions.

[Albert Chan]

(07-28-2019 12:02 AM)Albert Chan Wrote:  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

Note: above Σk^6 formula is from 0 to n = B(7,n)/7 + n^6

We may use symmetry, and do 0^6 to 3^6, then flip it, to get values (-3)^6 .. (+3)^6
This keep table numbers smaller, thus smaller forward difference numbers:

\(x^6 =
729\binom{x+3}{0}
-665\binom{x+3}{1}
+602\binom{x+3}{2}
-540\binom{x+3}{3}
+480\binom{x+3}{4}
-360\binom{x+3}{5}
+720\binom{x+3}{6}\)

Integrate (actually sum, but look like integration)

\(\large \frac{B^7}{7} \normalsize =
729\binom{x+3}{1}
-665\binom{x+3}{2}
+602\binom{x+3}{3}
-540\binom{x+3}{4}
+480\binom{x+3}{5}
-360\binom{x+3}{6}
+720\binom{x+3}{7} \) + C

Differentiate, evaluated at x=0, we get B(6) = RHS constant term.

First half of terms is not as easy to evaluate, because there was no factor x
(with a factor x, derivative at x=0 make the term goes away!)

But, it may still worth it, to keep forward difference table entries small.
Note: for B(m), n = floor(m/2), k = 1 .. n

\(\displaystyle \binom{x+n}{k} =
\left(\frac{x+n}{k}\right) \binom{x+n-1}{k-1}
\)

\(\displaystyle \binom{x+n}{k}' \bigg\rvert_{x=0} =
\left(\frac{1}{k}\right) \binom{n-1}{k-1} +
\left(\frac{n}{k}\right) \binom{x+n-1}{k-1}'\bigg\rvert_{x=0}
\)

I had managed to undo recursion, into iterative sums (written prove welcome!)

\(\displaystyle \binom{x+n}{k}' \bigg\rvert_{x=0}
\;=\; \sum_{j=1}^{k} \frac{(-1)^{j+1}}{j}\; \binom{n}{k-j}
\;=\; \binom{n}{k}\;\sum_{j=n-(k-1)}^{n} \frac{1}{j}
\)

Note that coefficient formula good for half terms (the other half is easy, see code)

Code:

from gmpy2 import mpq, mpz, xmpz
from itertools import izip, count

def bernoulli(m):
    if m == 0: return 1
    lst = [0, 1] + [ mpz(i)**m for i in range(2, (m+3)//2) ]
    lst = [(-x if m&1 else x) for x in lst[-1-(m&1):0:-1]] + lst
    lst = map(xmpz, lst)        # speedup fdiff_inplace
    fdiff_inplace(lst)          # forward diff table
    return sum(c*d for (c,d) in izip(bernoulli_coef(m//2), lst))

def fdiff_inplace(f):
    n = len(f) - 1
    for i in xrange(n):         # n forward diff(s)
        for j in xrange(n,i,-1):
            f[j] -= f[j-1]

def bernoulli_coef(n, mpq=mpq):
    'diff(comb(x+n,k)), k=1,2,3,..., evaluate at x=0'
    c, s, N = mpz(1), 0, n+1
    for j in xrange(n, 0, -1):  # tricky half
        c = c * j // (N-j)
        s += mpq(1, j)
        yield c * s
    s = mpq(1, N)
    for c in count(1):          # easy half
        yield s
        s *= mpq(-c, N+c)

I did not optimize away B(odd m), just to see if it matches ... it does! B(1) = -1/2, B(odd>1) = 0

>>> for m in range(2,31,2): print m, bernoulli(m)
...
2      1/6
4      -1/30
6      1/42
8      -1/30
10     5/66
12     -691/2730
14     7/6
16     -3617/510
18     43867/798
20     -174611/330
22     854513/138
24     -236364091/2730
26     8553103/6
28     -23749461029L/870
30     8615841276005L/14322

[Albert Chan]

Bernoulli number related links ...

Compare with forward difference table, Stirling numbers 2nd kind numbers are smaller
see thread, (HP15C)(HP67)(HP41C) Bernoulli Polynomials

Euler zigzag numbers (A) smaller still.
B(m) = A(m-1) * (-1)^(m/2-1) * m / (2^m*(2^m-1))
B(6) = 16 * (-1)^2 * 6 / (2^6*(2^6-1)) = 1/42

see thread, (28 48 49 50) Bernoulli Numbers

BCMATH programs for some number theory functions

CalcBn V3.0 from The Bernoulli Number Page
This use zeta function correction, same as this thread, but *much* faster!
On my laptop, it get B(100,000) exact fraction in 12 second!