(12C+) Bernoulli Number
07-27-2019, 06:41 AM (This post was last modified: 07-28-2019 06:08 AM by Gamo.)
Post: #1
 Gamo Senior Member Posts: 756 Joined: Dec 2016
(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

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, 12:41 PM (This post was last modified: 07-27-2019 10:41 PM by Albert Chan.)
Post: #2
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
(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)
07-27-2019, 01:40 PM
Post: #3
 Gamo Senior Member Posts: 756 Joined: Dec 2016
RE: (12C+) Bernoulli Number
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
07-27-2019, 07:49 PM
Post: #4
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (12C+) Bernoulli Number
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.
07-28-2019, 12:02 AM (This post was last modified: 08-25-2019 02:59 PM by Albert Chan.)
Post: #5
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
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

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
07-28-2019, 01:08 AM (This post was last modified: 07-28-2019 05:42 PM by Albert Chan.)
Post: #6
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
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
07-28-2019, 02:29 AM (This post was last modified: 07-28-2019 06:40 AM by Gamo.)
Post: #7
 Gamo Senior Member Posts: 756 Joined: Dec 2016
RE: (12C+) Bernoulli Number
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, 11:21 AM (This post was last modified: 07-18-2023 12:09 AM by John Keith.)
Post: #8
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (12C+) Bernoulli Number
(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

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.
07-31-2019, 05:14 PM (This post was last modified: 07-31-2019 11:04 PM by Albert Chan.)
Post: #9
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
Knowing the pattern of sk(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), s6(n) = n^7/7 - n^6/2 + n^5/2 + a*n^3 + b*n, where a, b are unknown

→ a*n^2 + b = s6(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
08-30-2023, 09:46 PM (This post was last modified: 08-30-2023 10:58 PM by Albert Chan.)
Post: #10
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
(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.
09-11-2023, 03:48 PM (This post was last modified: 09-12-2023 05:49 PM by Albert Chan.)
Post: #11
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
(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
09-12-2023, 05:59 PM
Post: #12
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (12C+) Bernoulli Number
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!
 « Next Oldest | Next Newest »

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