+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (HP15C)(HP67)(HP41C) Bernoulli Polynomials (/thread-20416.html)
Pages: 1 2
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Paul Dale - 09-03-2023 05:35 AM
How does the WP 34S fare for these?
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 09-05-2023 05:57 PM
(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.
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-05-2023 07:09 PM
(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?
![[Image: 463600a23eacbff5cf37ab115762123b6186d5bc]](https://wikimedia.org/api/rest_v1/media/math/render/svg/463600a23eacbff5cf37ab115762123b6186d5bc)
The algorithm is simply rewriting B(2n) as function of zeta(2n),
and B(2n) denominator = product(p, such that (p-1)|2n)
![[Image: e5fed8a1d28f2196fc17d7631657be9d35cedb85]](https://wikimedia.org/api/rest_v1/media/math/render/svg/e5fed8a1d28f2196fc17d7631657be9d35cedb85)
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, ...
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-05-2023 08:15 PM
(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
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 09-05-2023 09:47 PM
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.
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - Albert Chan - 09-05-2023 10:39 PM
(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
endlua> 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.
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 09-07-2023 04:31 PM
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.
RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials - John Keith - 09-08-2023 11:26 PM
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
\>>
\>>