(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[Albert Chan]

The reason we sidetrack for Bernoulli numbers is we wanted to use this:



With previous post for bernoulli numbers (O(1), speed and storage), we have:

Code:

function b0(m)
    if m%2==1 then return m==1 and -1/2 or 0 end
    if m == 0 then return 1 end
    local b = (-1)^(m/2+1) * 2*tgamma(m+1)/pi^m
    local d = 2^m-1
    local n = floor(b + b*3^-m)
    return (n+0.5)/d, b/(d+1) -- ratio = zeta(m)
end

Code:

function b(m, x)
    local c, t = 1, 1
    for k = 1, m do
        c = c*(m-k+1)/k
        t = t*x + c*b0(k)
    end
    return t
end

lua> b(3, 5.5)
123.75
lua> b(5, 3.5)
220.9375

lua> b(15, pi)
640433.1973240786
lua> b(16, pi)
1462871.9320025162

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!)
Also, I remove complex number math for sign of b0(m)

Update 3: return bad b0(m) too, before zeta correction. Ratio is zeta(m)

lua> good, bad = b0(2)
lua> good, bad, good/bad -- = zeta(2)
0.16666666666666666      0.10132118364233778      1.6449340668482262