(28 48 49 50) Bernoulli Numbers
|
09-05-2023, 05:27 PM
(This post was last modified: 09-10-2023 07:55 PM by John Keith.)
Post: #1
|
|||
|
|||
(28 48 49 50) Bernoulli Numbers
Note:The programs in this post return Bernoulli numbers as approximate (real) numbers. In post #6 below is a program for the HP-28 and HP-48 that returns the exact numerator and denominator for all m <= 28. Post #12 below has Exact mode programs for the HP 49/50 which are significantly faster than IBERNOULLI.
Given an integer n on the stack, this program returns the nth Bernoulli number as a real number. The result is accurate to 12 digits for numbers up to B(18) and to 11 or 12 digits up to B(40). The largest error is 5 ULP's for B(20), all others are 2 ULP's or less. The program has not been tested for n > 40 but larger Bernoulli numbers are > 10^16 and thus not very useful in most cases. This program grew out of the discussion in this thread and is based on this section of the Wikipedia page. Code:
This next program returns a list of Bernoulli numbers from 0 through n. The value of n is constrained to be an even number >= 4 in order to avoid making a large program even larger. For the HP-28, remove the NEWOB in line 9. Code:
|
|||
09-06-2023, 08:48 AM
Post: #2
|
|||
|
|||
RE: (48G) Bernoulli numbers
On a 50g this programme
Code: Size: 277. for input 1 100 returns ' -94598037819122125295227433069493721872702841533066936133385696204311395415197247711 /33330' in 42.2994s |
|||
09-06-2023, 11:04 AM
Post: #3
|
|||
|
|||
RE: (48G) Bernoulli numbers
Wow, that's fast! The built-in command takes about 77 seconds on my 50g. What algorithm are you using?
|
|||
09-06-2023, 02:34 PM
(This post was last modified: 09-06-2023 02:44 PM by Gerald H.)
Post: #4
|
|||
|
|||
RE: (48G) Bernoulli numbers
I suppose I wrote the programme when the 49G was introduced as built in
FPTR2 ^IBERNOULLI was so slow. I have no records from that time & can only guess I copied some well known algorithm. The programme actually gives the polynomial for input of 'x' 100 if that helps (shouldn't try, takes too long). |
|||
09-07-2023, 03:37 PM
Post: #5
|
|||
|
|||
RE: (48G) Bernoulli numbers
Here is lua function for zigzag numbers, build 1 zigzag at a time
Code: function zigzag(t) -- Euler zigzag numbers lua> m = 16 lua> L = m/2-1 -- 1(excluded) to m-1 step 2, zigzags to reach = ((m-1)-1)/2 = m/2-1 lua> t = {[0] = 1} -- t[0], t[-0] = T(0), T(1) lua> for k=1,L do t=zigzag(t) end lua> t[-#t] -- = T(m-1) 1903757312 lua> p = 2^m lua> n = (-1)^L * m * t[-#t] * 2/p lua> d = 2*(p-1) -- B(m) denominator = 2*odd lua> n, d, n/d -929569 131070 -7.092156862745098 Array t can be extended, since biggest element, T(m-1), does not overflow yet. lua> m = 100 lua> t[-#t] = require'mapm'.new(t[-#t]) -- pollute array with mapm number lua> L2 = m/2-1 lua> for k=1,L2-L do t=zigzag(t) end lua> t[-#t] -- = T(m-1) 4.560851661680111182104382953145169718558E+136 -- Note: mapm +/-/* is exact, above holds *exact* T(m-1) -- mapm.gcd use efficient binary gcd algorithm. I don't bother with denominator = 2*odd lua> ;mapm.places(136) -- T(m-1) exponent more than enough precision lua> n = (-1)^L2 * m * t[-#t] lua> d = mapm.pow(2, m) lua> d = d*(d-1) lua> g = mapm.gcd(n, d) lua> n, d = n/g, d/g lua> n:tofixed(0) .. ' / ' .. d:tofixed(0) -- B(m) reduced fraction -94598037819122125295227433069493721872702841533066936133385696204311395415197247711 / 33330 |
|||
09-07-2023, 04:18 PM
Post: #6
|
|||
|
|||
RE: (48G) Bernoulli numbers
The following program is a translation of Albert Chan's Lua program here. Given m on the stack, the program returns the numerator on level 2 and the denominator on level 1. The results are exact for all m <= 28. The program is for the HP-28 and HP-48. It can be used in approximate mode on the HP49 and 50 but is not really necessary because those models have IBERNOULLI built in, and can use Gerald's faster program above.
The program returns 1, 1 for B(0) and 0, 1 for odd m > 2 so that the program will always return exactly two numbers regardless of input. The numerators will not be correct for n > 28 but the denominators will be. Code:
|
|||
09-08-2023, 08:09 PM
Post: #7
|
|||
|
|||
RE: (28 48) Bernoulli numbers
Update: The programs in the first post have been replaced by new ones that are faster as well as HP-28 and HP-48S compatible. The thread title has been changed to reflect these changes. The program in post 6 above is still preferable in most cases as it returns the Bernoulli number <= B(28) as a fraction with separate numerator and denominator.
|
|||
09-09-2023, 06:33 PM
Post: #8
|
|||
|
|||
RE: (28 48) Bernoulli numbers
(09-07-2023 04:18 PM)John Keith Wrote: The following program is a translation of Albert Chan's Lua program here. I have reduced getting reduced B(even m) denominator, from O(m) to O(√m) Let's show with example, m = 30 Let m = (p-1)*(q-1), where bolded (p, q) are primes. divisors(m) = 1, 2, 3, 5, 6, 10, 15, 30 divisors(m)+1 = p = 2, 3, 4, 6, 7, 11, 16, 31 m/(p-1) + 1 = q = 31, 16, 11, 7, 6, 4, 3, 2 We only need half the p's, and corresponding q's is p's other half. d = product(prime p, (p-1) | m) = (half p's) * (half q's) = (2*3) * (7*11*31) = 14322 By design (blue region), p factors <= q factors We keep 1 factor, if p=q, to avoid double counting. Since we start off half p's first, half q's loop will quit if p >= q. Code: function min_d(m) -- B(even m>=2) reduced denominator lua> for m=2, 40, 2 do print(m, min_d(m)) end 2 6 4 30 6 42 8 30 10 66 12 2730 14 6 16 510 18 798 20 330 22 138 24 2730 26 6 28 870 30 14322 32 510 34 6 36 1919190 38 6 40 13530 |
|||
09-09-2023, 07:34 PM
(This post was last modified: 09-10-2023 05:33 PM by Albert Chan.)
Post: #9
|
|||
|
|||
RE: (28 48) Bernoulli numbers
(09-09-2023 06:33 PM)Albert Chan Wrote: By design (blue region), p factors <= q factors Simpler version, removed test for p >= q We know halfp loops, n = 2^m-1 had halfp factors completely removed. Test for prime if q=p, n%q==0 is false, thus no double counting factors. max(p factors) ≤ min(q factors) hi + 1 ≤ m / hi + 1 hi^2 ≤ m hi ≤ √m Code: function min_d(m) -- B(even m>=2) reduced denominator if we set hi = m, and remove k loops, we get back O(m) version. Benchmarks suggested O(√m) version is faster for even m ≥ 12 |
|||
09-10-2023, 03:24 PM
Post: #10
|
|||
|
|||
RE: (28 48) Bernoulli numbers
That's very neat but unfortunately it turns out to be slower than your original, smaller program for the limited range of Bernoulli numbers that can be computed this way on the HP-28/48. I imagine that it would be an improvement for systems with multiple precision floats where m can be much larger.
|
|||
09-10-2023, 07:39 PM
(This post was last modified: 09-10-2023 11:17 PM by Albert Chan.)
Post: #11
|
|||
|
|||
RE: (28 48) Bernoulli numbers
Hi, John Keith
I was also surprised that O(m) to O(√m) does not translate to much speed gain. Yes, O(√m) version eventually is faster, but not by much. The reason is D = 2*(2^m-1) factors between 2 and m+1 matched d factors, with few false positives. This make O(m) code very efficient. (it skipped D nonfactors) m = 2 to 52 step 2, there are only 3 m's that have false positives. (bold) m = 24: d = 2*3*5*7*13 D = 2*3²*5*7*13*17*241 m = 40: d = 2*3*5*11*41 D = 2*3*5²*11*17*31*41*61681 m = 50: d = 2*3*11 D = 2*3*11*31*251*601*1801*4051 Instead of half q's, we may use a few q's, to reduce search space. With 1 q, it already beat O(√m) version, all the way to m=52 (IEEE double limit) Code: function min_d(m) -- B(even m>=2) reduced denominator |
|||
09-10-2023, 07:45 PM
(This post was last modified: 09-10-2023 07:49 PM by John Keith.)
Post: #12
|
|||
|
|||
RE: (28 48) Bernoulli numbers
Albert Chan's program (see post # 6 above) is incredibly fast for numbers up to B(28) but can't be used for larger numbers, even on the HP 49 and 50, because it requires floating-point numbers. However, it turns out that the method using tangent numbers ( see post # 1) is significantly faster than the IBERNOULLI function in the 49/50, and with exact integers can compute Bernoulli numbers as large as memory allows.
The following program, similar to the first program in post #1, computes the nth Bernoulli number. It is over 3 * as fast as IBERNOULLI for B(100). The program requires GoferLists and must be used in Exact mode. Code:
The next program, similar to the second program in post #1, returns a list of Bernoulli numbers from 0..n. It is about 8 * as fast as a loop using IBERNOULLI. Exact mode required. Code:
Additional program notes: Though ListExt is newer, more flexible and generally faster than GoferLists, I am using GoferLists here because the commands Scanr and Last make the programs shorter and more readable. No special code is necessary to determine the denominator because exact integer division automatically reduces fractions to lowest terms. |
|||
09-17-2023, 02:32 PM
Post: #13
|
|||
|
|||
RE: (28 48 49 50) Bernoulli Numbers
On a 50g this revised programme
Code: Size: 251. for input 102 returns '3204019410860907078243020782116241775491817197152717450679002501086861530836678158791/4326' in 25.1814s |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 5 Guest(s)