A fast Bernoulli Number method for the HP Prime
Message #1 Posted by Namir on 21 Nov 2013, 9:48 a.m.
During HHC2013 I mentioned at the end of my presentation that the HP Prime had a wonderful collection of functions but lacked functions like the Bessel functions, Bernoulli numbers, and Bernoulli polynomial. I commented that calculating the latter two was not easy, and Cyril affirmed my observation.
Well I have been tinkering with calculating the Bernoulli numbers and Bernoulli polynomials. The Bernoulli numbers have the following features:
1. Aside from Bern(1), all Bernoulli numbers for odd values are zero.
2. The Bernoulli numbers for even values are not zero and alternate in signs.
3. Bernoulli numbers start as having small values that increase to very large values. I suspect that this variation in values causes stability problems in calculating Bernoulli numbers using certain algorithms.
Starting with the wealth of equations for the Bernoulli numbers in Wikipedia, I began trying various equations presented there (and I recommend you look at Wikipedia's web pages). After implementing various algorithms, observing the stability of their results, I found the following:
1. The recursive method were dead on in giving accurate results.
2. The series approximation method and the two matrix-based equations gave very good accuracy--not perfect but very good for the computation effort and code involved.
3. The rest of methods, like double summations, the Akiyama–Tanigawa method, and a variant of the Akiyama–Tanigawa method that I stumbled upon on the Internet, worked well for up to Bernoulli(10) and then gave bad results for higher numbers!!!
I implemented one of the matrix methods that give results that are close to the exact values. The code enjoys simplicity compared to implementing other algorithms (which I did in Excel VBA):
EXPORT BernMat(n)
BEGIN
LOCAL i, j, m, fact, mat;
// NOTE: Replace the character # with the proper
// not-equal character for the code to run correctly
IF n#1 AND (n MOD 2) == 1 THEN
RETURN 0;
END;
IF n== 1 THEN
RETURN -0.5;
END;
IF n==2 THEN
RETURN 1.0/6.0;
END;
m:=n+1;
// Create the n+1 by n+1 matrix
mat:=MAKEMAT(0,m,m);
// only non-zero value in column n+1
mat(1, m):=1;
// Calculate elements of the first matrix column
// Fact store the factorials of 2 to n+1
fact:=1;
FOR i from 1 TO m DO
fact:=fact*i;
mat(i, 1):=1/fact;
END;
// copy matrix elements to fill columns 2 to n
FOR j FROM 2 TO n DO
FOR i FROM j TO m DO
mat(i, j):=mat(i-1, j-1);
END;
END;
// divide Fact by m to obtain n!
RETURN fact / m * DET(mat);
END;
Calculating the coefficients of the Bernoulli polynomial is very simple once you have the Bernoulli numbers!!
Enjoy!!
Namir
Edited: 21 Nov 2013, 3:23 p.m. after one or more responses were posted
| 
Go back to the main exhibit hall