(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[Albert Chan]

(09-05-2023 09:47 PM)John Keith Wrote:  Of course, we then need a way to compute divisors of m and primes or GCD, which no HP calculator before the HP49 has. There are good factoring libraries for the HP-48 but this would not be easy for keystroke programmable calculators.

We can use the fact that D = 2*(2^m-1) has all d factors.
Also, d prime factors ≤ (m+1), relatively small number.

Code:

function ratioB(m)
    local b = I.real(-2*tgamma(m+1)/(2*pi*I)^m)
    local n, d = 2^m-1, 1
    for k=3, m+1, 2 do 
        if n%k ~= 0 then continue end
        repeat n=n/k until n%k ~= 0
        if m%(k-1)==0 then d = d*k end
    end
    n = floor(b*d * (1+2^-m+3^-m))
    return 2*n+1, 2*d    
end

lua> for m=2,40,2 do print(m, ratioB(m)) end

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  -23749461029            870
30  8615841276005           14322
32  -7709321041217          510
34  2577687858367           6
36  -26315271553053516000   1919190
38  2929993913841563        6
40  -261082718496449660000  13530

Note that all B(m) denominator's are correct.
IEEE double (53 bits) managed to get upto B(34) fraction exactly.