(HP15C)(HP67)(HP41C) Bernoulli Polynomials

09032023, 05:35 AM
Post: #21




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
How does the WP 34S fare for these?


09052023, 05:57 PM
Post: #22




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(09022023 05:44 PM)Albert Chan Wrote: We may be able to correct for tiny errors. B(m) largest denominator is 2*(2^m1). 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 HP48 (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. 

09052023, 07:09 PM
(This post was last modified: 09112023 08:13 PM by Albert Chan.)
Post: #23




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(09052023 05:57 PM)John Keith Wrote: This is quite interesting, I haven't seen this method before. The algorithm is simply rewriting B(2n) as function of zeta(2n), and B(2n) denominator = product(p, such that (p1)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 (09032023 01:20 AM)Albert Chan Wrote: Update 1: halved d, and use floor instead of round 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^m1 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^m1 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, ... 

09052023, 08:15 PM
Post: #24




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(09022023 05:44 PM)Albert Chan Wrote: Denominator likely smaller, but required work. d = product(p, for (p1)  even m) Proof: B(even m) denominator d, divides D = 2*(2^m1) If p=2, (p1)=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 (p1)  m, then m = k*(p1) 2^m  1 ≡ (2^(p1))^k  1 ≡ 1^k  1 ≡ 0 (mod p) // Fermat's little theorem D have all d prime factors > d divides D 

09052023, 09:47 PM
(This post was last modified: 09052023 09:48 PM by John Keith.)
Post: #25




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
Using the formula d = Product p, (p1)  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 HP48 but this would not be easy for keystroke programmable calculators.


09052023, 10:39 PM
Post: #26




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
(09052023 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 HP48 but this would not be easy for keystroke programmable calculators. We can use the fact that D = 2*(2^m1) has all d factors. Also, d prime factors ≤ (m+1), relatively small number. Code: function ratioB(m) 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. 

09072023, 04:31 PM
Post: #27




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
That's amazing, Albert! The program gives exact results on 12digit calculators up to B(28). I have posted a translation for the HP28 and HP48 here.
This program should be excellent for the DM42 and Free42 which have 34digit numbers, and even better on NewRPL which has variable precision up to 2000 digits. 

09082023, 11:26 PM
Post: #28




RE: (HP15C)(HP67)(HP41C) Bernoulli Polynomials
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:


« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)