HP Forums
A fast Bernoulli Number method for the HP Prime - Printable Version

+- HP Forums (https://archived.hpcalc.org/museumforum)
+-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html)
+--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html)
+--- Thread: A fast Bernoulli Number method for the HP Prime (/thread-256591.html)



A fast Bernoulli Number method for the HP Prime - Namir - 11-21-2013

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


Re: A fast Bernoulli Number method for the HP Prime - Alberto Candel - 11-21-2013

Hi,
If all that you need is an approximation to B(n), have you considered using Euler's formula relating Zeta(2n) to B(2n)? The Prime has a Zeta and computes its values at the even integers rather fast. They come out as a rational number times a power of pi. I do not know what algorithm is used for this computation.


Re: A fast Bernoulli Number method for the HP Prime - parisse - 11-21-2013

It computes the Bernoulli number, the command bernoulli is available in Xcas but for some obscure reason it is not on Prime. The code uses the recurrence relation, it stores intermediate result in a list to avoid recursive calls:

gen bernoulli(const gen & x){
if ( (x.type!=_INT_) || (x.val<0))// Should add bernoulli polynomials
return gensizeerr(gettext("bernoulli"));
int n=x.val;
if (!n)
return plus_one;
if (n==1)
return minus_one_half;
if (n%2)
return zero;
gen a(plus_one);
gen b(rdiv(1-n,plus_two,context0));
vecteur bi(makevecteur(plus_one,minus_one_half));
int i=2;
for (; i< n-1; i+=2){
// compute bernoulli(i)
gen A=1;
gen B=gen(1-i)/2;
for (int j=2; j<i-1;j+=2){
A=iquo( A*gen(i+3-j)*gen(i+2-j),(j-1)*j);
B=B+A* bi[j];
}
bi.push_back(-B/gen(i+1));
bi.push_back(0);
a=iquo( (a*gen(n+3-i)*gen(n+2-i)),((i-1)*i));
b=b+a* bi[i];
}
return rdiv(-b,n+1,context0);
}


Re: A fast Bernoulli Number method for the HP Prime - Han - 11-21-2013

Quote:
  gen bernoulli(const gen & x){
if ( (x.type!=_INT_) || (x.val<0))// Should add bernoulli polynomials
return gensizeerr(gettext("bernoulli"));
int n=x.val;
if (!n)
return plus_one;
if (n==1)
return minus_one_half;
if (n%2)
return zero;
gen a(plus_one);
gen b(rdiv(1-n,plus_two,context0));
vecteur bi(makevecteur(plus_one,minus_one_half));
int i=2;
for (; i< n-1; i+=2){
// compute bernoulli(i)
gen A=1;
gen B=gen(1-i)/2;
for (int j=2; j<i-1;j+=2){
A=iquo( A*gen(i+3-j)*gen(i+2-j),(j-1)*j);
B=B+A* bi[j];
}
bi.push_back(-B/gen(i+1));
bi.push_back(0);
a=iquo( (a*gen(n+3-i)*gen(n+2-i)),((i-1)*i));
b=b+a* bi[i];
}
return rdiv(-b,n+1,context0);
}

Added formatting to help make it more legible.


Re: A fast Bernoulli Number method for the HP Prime - Ángel Martin - 11-21-2013

That's also how I did it on the SandMath - see below.

n * Zeta (1-n) = - Bn

01	LBL "BN2" 
02 X=1?
03 GTO 01
04 X=0?
05 INCX
06 X=1?
07 RTN
08 ODD?
09 CLX
10 X=0?
11 RTN
12 2
13 X<Y?
14 GTO 00
15 6
16 1/X
17 RTN
18 LBL 00
19 ST+ X
20 X<Y?
21 GTO 00
22 -30
23 1/X
24 RTN
25 LBL 00
26 X<>Y
27 LBL 01
28 STO M
29 CHS
30 INCX
31 ZETA
32 RCL M
33 CHS
34 *
35 END

Edited: 21 Nov 2013, 11:48 a.m.


Re: A fast Bernoulli Number method for the HP Prime - Alberto Candel - 11-21-2013

If I understand correctly, Zeta(2n) is computed via Euler's identity Zeta(2n)= (+-)B(2n)*(2pi)^2n /2(2n!). How are other values of Zeta computed? Like Apery's Zeta(3). Thanks


Re: A fast Bernoulli Number method for the HP Prime - Namir - 11-21-2013

Alberto,

Check Wikipedia's Bernoulli Number web page and look for the expressions for Bn that use matrices.

Namir


Re: A fast Bernoulli Number method for the HP Prime - Paul Dale - 11-21-2013

The 34S uses its implementation of the Zeta function to calculate the Bernoulli numbers.


- Pauli


Re: A fast Bernoulli Number method for the HP Prime - Namir - 11-21-2013

Here is the HP Prime program that calculates Bernoulli number using series approximation:

EXPORT BernSer(n,toler)
BEGIN
LOCAL sum,term,i,chs;

IF n==1 THEN
RETURN -0.5;
END;

IF n Mod 2 == 1 THEN
RETURN 0;
END;

sum:=1;
i:=2;
term:=2*toler;
WHILE term>toler DO
term:=1/i^n;
sum:=sum+term;
i:=i+1;
END;

IF n MOD 4 == 0 THEN
chs:=-1;
ELSE
chs=1;
END;
RETURN chs*(n!)/pi^n/2^(n-1)*sum;
END;

NOTE: I updated the listing tweak it.


Edited: 23 Nov 2013, 10:41 a.m. after one or more responses were posted


Re: A fast Bernoulli Number method for the HP Prime - Joe Horn - 11-22-2013

Quote:
That's also how I did it on the SandMath - see below.

n * Zeta (1-n) = - Bn


Cool. So, on Prime, define this in CAS:

Bernoulli(x) := –x*Zeta(1–x)

Cranks out Bernoulli(200), all 224 digits over 7 digits, in about 3.7 seconds on the hardware.

Edit: The above gets the wrong answer for Bernoulli(0). Gotta add that as an exception. It should return 1.

Edit 2: Bernoulli(x) := IFTE(x=0,1,–x*Zeta(1–x))

Edited: 22 Nov 2013, 7:55 a.m.


Re: A fast Bernoulli Number method for the HP Prime - Gerson W. Barbosa - 11-22-2013

Great! This works also on my HP 50g with this Riemann's Zeta function implementation:

<< -199 Zeta 200 * NEG >> TEVAL     -->  -3.64707726448E215   (  0.47 s)

<< 200 IBERNOULLI -> NUM >> TEVAL --> -3.64707726452E215 (486.83 s)
The HP-28S version returns -3.64707724956E215 in a couple of seconds.




Re: A fast Bernoulli Number method for the HP Prime - Namir - 11-22-2013

Joe,

The availability of the Zeta function makes evaluating Bernoulli numbers and polynomials very trivial. The code I presented on this site calculates the Bernoulli numbers using simpler functions (well .. if you can call the matrix determinant simple ... :-)). My second listing uses more basic operations to calculate the Bernoulli numbers. I should be able to post a listing for the HP-41C soon (I am sure J-m Balliard has done that already)!! Even an HP-67 version is doable.

It seems that the HP Prime and the wp34s calculators share this advantage!!
Namir


Re: A fast Bernoulli Number method for the HP Prime - Alberto Candel - 11-22-2013

Hi all,
Does anyone have a reference to how the HP's (prime, 50, 34s) compute Zeta(s)? Because, paraphrasing Namir, using zeta to compute B(n) is kind of a blind roundabout.

-a

[added] my bad. I just saw some of the code for zeta in the post by Gerson above. Thanks
[/added]

Edited: 22 Nov 2013, 11:20 a.m.


Re: A fast Bernoulli Number method for the HP Prime - Gerson W. Barbosa - 11-22-2013

Hi Alberto,

If you've followed the links you may have found another method, Borwein's algorithm as implemented by Thomas Ritschel in message #5 here:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv021.cgi?read=235199

That one is better, but I chose the older method because it was easier to follow in the reference I had (message #20 in this thread). I've used fix parameters that grant accuracy in the real domain only. Arguments with absolute imaginary parts greater than 50 or so would require more terms. That was enough for my needs, but it's not a definitive zeta implementation for complex arguments yet.

Gerson.


Edited: 22 Nov 2013, 1:41 p.m.


Re: A fast Bernoulli Number method for the HP Prime - Namir - 11-22-2013

And here is the HP-41C implementation for Bernoulli numbers using series approximation. It requires no modules or synthetic commands:

LBL "BERSER"
LBL A
CF 22
"TOLER?"
PROMPT
FS? 22
STO 01 # store new tolerance value if user entered it
"N?"
PROMPT
STO 00
1
X<>Y
X=Y? # n=1?
GTO 11
2
MOD
1
X=Y? # is n odd?
GTO 10
1
STO 03 # sum = 1
STO 02 # i = 1
LBL 01 # Start summation loop
1
STO+ 02 # i=i+1
RCL 02
RCL 00
Y^X
1/X # Calculate 1/i^n
STO+ 03 # sum = sum + 1/i^n
RCL 00
X<=Y? # toler <= term?
GTO 01 # resume to the next iteration
CF 00 # start calculating other factor
RCL 00
4
MOD
X=0?
SF 00
RCL 00
FACT # n!
PI
RCL 00
Y^X # pi^n
/
2
RCL 00
1
-
Y^X # 2^(n-1)
/
RCL 03
*
FS?C 00 # change sign?
CHS
RTN
LBL 10
0
RTN
LBL 11
-0.5
RTN

Enjoy!


Edited: 23 Nov 2013, 10:52 a.m.


Re: A fast Bernoulli Number method for the HP Prime - Alberto Candel - 11-22-2013

Hi Namir,

Thanks, nice, just tried it (I don't have a 41 thou). Is there a way to modify the code so as to output the actual number as a fraction?

-a


Re: HP Prime numerical restrictions? - Alasdair McAndrew - 11-25-2013

There are no known closed form expressions for odd values of the zeta function. Apery was able to prove that zeta(3) was irrational, and that's pretty much as far as people know about zeta(2n+1). According to the wikipedia page, infinite values of zeta(2n+1) are irrational, and at least one of zeta(5), zeta(7), zeta(9), zeta(11) is irrational. But no closed forms!