HP Forums - (12C+) Bernoulli Number

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

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

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

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.

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 Smile

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

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

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

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

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

(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

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!