(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[Albert Chan]

(09-05-2023 05:57 PM)John Keith Wrote:  This is quite interesting, I haven't seen this method before.
Can it be modified to avoid complex numbers so that it can be used on older and simpler calculators?

The algorithm is simply rewriting B(2n) as function of zeta(2n),
and B(2n) denominator = product(p, such that (p-1)|2n)



It does not involve complex numbers. Complex number form is just more compact.

B(2n) = Re(-2 (2n)! / (2*pi*i)^(2n)) * zeta(2n)

(2*pi*i)^(2n) = (2*pi)^(2n) * i^(2n) = (2*pi)^(2n) * (-1)^n

(09-03-2023 01:20 AM)Albert Chan Wrote:  Update 1: halved d, and use floor instead of round

b0(m) = (n+0.5)/d = (2n+1)/(2d) = odd/even, as expected

Update 2: zeta = (1+3^-m)/(1-2^-m) estimate is better than (1+2^-m+3^-m)
We could implement this without cost (actually, cheaper!)

Redo the same examples, with update #1, #2 ideas.
Here, we add sign afterwards, but it could be added to b beforehand.

lua> m = 12
lua> b, d = 2*fac(m)/pi^m, 2^m-1
lua> b, d
1036.4980453216303      4095
lua> _ + b*3^-m
1036.4999956755648
lua> - 1036.5 / 4095 -- = B(12)
-0.2531135531135531

lua> m = 16
lua> b, d = 2*fac(m)/pi^m, 2^m-1
lua> b, d
464784.48919973      65535
lua> _ + b*3^-m
464784.49999694
lua> - 464784.5 / 65535 -- = B(16)
-7.092156862745098

Note: next corrections (if needed) are b*5^-m, b*7^-m, ...