(HP15C)(HP67)(HP41C) Bernoulli Polynomials

[Paul Dale]

How does the WP 34S fare for these?

[John Keith]

(09-02-2023 05:44 PM)Albert Chan Wrote:  We may be able to correct for tiny errors. B(m) largest denominator is 2*(2^m-1).

d = 2*(2^12-1) = 8190 (largest d)
d * -0.25312 = -2073.0528            → B(12) = -2073/8190 = -691/2730

Denominator likely smaller, but required work. d = product(p, for (p-1) | even m)
This is why d is divisible by 6. 2-1=1, 3-1=2, both divides even m.

d = 2*3*5*7*13 = 2730
d * -0.25312 = -691.0176             → B(12) = -691/2730

---

Amazingly, we get correct B(m = 12), without correction (algorithm based on above link)
If m divisible by 4, B(m) is negative. (positive otherwise)

lua> m = 12
lua> d = 2*(2^m-1)
lua> b = 4*d/(d+2) * fac(m)/pi^m
lua> 1+b, d
2073.4899880820685      8190
lua> _ + b*2^-m
2073.995967083065

--> B(12) = -2073/8190 = -691/2730

lua> m = 16
lua> d = 2*(2^m-1)
lua> b = 4*d/(d+2) * fac(m)/pi^m
lua> 1+b, d
929555.7943024995      131070
lua> _ + b*2^-m
929569.9781830278
lua> _ + b*3^-m
929569.9997771184

--> B(16) = -929569/131070 = -3617/510

---

Rounding errors might cause off-by-1 error with B(m) numerator.
It may be safer to drop numbers by 1/2, equivalent to rounding without 1+

Code:

gcd = require'factor'.gcd
fac = fn'x: tgamma(x+1)'
roughB = fn'm: I.real(-2*fac(m)/(2*pi*I)^m) * (1+2^-m+3^-m)'
ratioB = fn'm,n,d: d=2*(2^m-1); n=round(roughB(m)*d); m=gcd(n,d); n/m, d/m'

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?

Additionally, I have tested the method that I referred to above using zigzag (actually tangent) numbers, and the accuracy is far better than the method using Stirling numbers and factorials. On the HP-48 (12 digits) the results are correct to all 12 digits up to B(18) and to 11 or 12 digits up to B(40), which is as far as I tested. Programs are here. Also see this paper for further information and algorithms.

[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, ...

[Albert Chan]

(09-02-2023 05:44 PM)Albert Chan Wrote:  Denominator likely smaller, but required work. d = product(p, for (p-1) | even m)

Proof: B(even m) denominator d, divides D = 2*(2^m-1)

If p=2, (p-1)=1, always divides m, and D also have factor 2.
With even denominator, B(even m) numerator must be odd.

All odd p, such that (p-1) | m, then m = k*(p-1)

2^m - 1 ≡ (2^(p-1))^k - 1 ≡ 1^k - 1 ≡ 0 (mod p)      // Fermat's little theorem

D have all d prime factors --> d divides D

[John Keith]

Using the formula d = Product p, (p-1) | m from your 2nd link in post #19 along with approximate values of B(n) computed by tangent numbers, we can get exact a up to B(28) with 12 digits, or B(26) with 10 digits. 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.

[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.

[John Keith]

That's amazing, Albert! The program gives exact results on 12-digit calculators up to B(28). I have posted a translation for the HP-28 and HP-48 here.

This program should be excellent for the DM42 and Free42 which have 34-digit numbers, and even better on NewRPL which has variable precision up to 2000 digits.

[John Keith]

Now that we have (mostly) conquered Bernoulli numbers, we can get back to the original purpose of this thread which is Bernoulli polynomials. Albert's second program in post #20 is simple and clearly written, and should be easy to implement on any HP calculator that can handle Bernoulli numbers.

To start off, here is an implementation that will run on any RPL calculator. In this example, BERN can refer to the first program in this post in the other thread, or the program here. In the latter case, BERN should be followed by / (divide).

Code:


\<< \-> m x
  \<< 1. 1. 1. m
    FOR k m k - 1. + * k /
      SWAP x * OVER k BERN * + SWAP
    NEXT DROP
  \>>
\>>