(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
 John Keith Senior Member Posts: 1,040 Joined: Dec 2013
(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:

\<< DUP 3.
IF <
THEN 1. -2. INV 6. INV 3. \->LIST          @ Initial values for n in [0..2]
SWAP 1. + GET
ELSE
IF DUP 2. MOD                            @ B(n) = 0 for odd n > 2
THEN DROP 0.
ELSE 2. / \-> n                          @ Even values only
\<< 0. 1. 2. n                         @ Loop to compute tangent number
FOR k 0. 3. k 2. *
FOR t DUP t ROLL +
NEXT 0. 3. k 2. * 1. +
FOR t DUP t ROLL +
NEXT
NEXT n 2. * ROLLD n 2. * 1. - DROPN  @ nth tangent number
n DUP 2. MOD :: NEG IFT *            @ Correct sign
4. n ^ DUP 1. - SWAP 2. / NEG * /    @ Divide by (2^n-4^n)/2
\>>
END
END
\>>

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:

\<< DUP 2. MOD - 4. MAX 2. / \-> n         @ Constrain to even n >= 4
\<< 1. -2. INV 6. INV { 0. 1. } 2. n     @ First 3 terms
FOR k 0. SWAP LIST\-> 2. + \-> s       @ Loop to compute tangent numbers
\<< 0. 3. s
FOR t DUP t ROLL +
NEXT 0. 3. s 1. +
FOR t DUP t ROLL +
NEXT s \->LIST
\>> DUP DUP SIZE GET NEWOB           @ Tangent number
k 2.
IF MOD
THEN NEG                             @ Correct sign
END k *                              @ multiply by k/2
4. k ^ DUP 1. - SWAP -2. / * / SWAP  @  Divide by (2^n-4^n)/2
NEXT DROP n 2. * 1. + \->LIST
\>>
\>>
09-06-2023, 08:48 AM
Post: #2
 Gerald H Senior Member Posts: 1,627 Joined: May 2014
RE: (48G) Bernoulli numbers
On a 50g this programme

Code:
Size: 277.

CkSum: # 1FA8h

::
CK2&Dispatch
# AFF
::
FPTR2 ^DupQIsZero?
case2drop
ZINT 1
FPTR2 ^DupZIsOne?
case
::
ZINT 2
FPTR2 ^NDXFext
FPTR2 ^QSub
FPTR2 ^CASCOMPEVAL
;
FPTR2 ^DupZIsEven?
NOT
3PICK
CRUNCH
%0=
AND
case2drop
ZINT 0
DUP
ZINT 2
3UNROLL
FPTR2 ^Z2BIN
ZINT 0
ZINT 0
ZINT -1
4PICK
#2+
ONE_DO
ZINT 1
7PICK
INDEX@
FPTR2 ^#>Z
FPTR2 ^QDiv
FPTR2 ^QSub
FPTR2 ^QMul
FPTR2 ^CASCOMPEVAL
5PICK
5PICK
FPTR2 ^PPow#
ROT
FPTR2 ^CASCOMPEVAL
SWAP
5ROLL
ZINT 1
5UNROLL
2DUP
FPTR2 ^QMul
INDEX@
FPTR2 ^#>Z
FPTR2 ^QDiv
4ROLL
FPTR2 ^CASCOMPEVAL
3UNROLL
LOOP
2DROP
4UNROLL3DROP
;
;

for input

1
100

returns

'
-94598037819122125295227433069493721872702841533066936133385696204311395415197247​711
/33330'

in

42.2994s
09-06-2023, 11:04 AM
Post: #3
 John Keith Senior Member Posts: 1,040 Joined: Dec 2013
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
 Gerald H Senior Member Posts: 1,627 Joined: May 2014
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
 Albert Chan Senior Member Posts: 2,705 Joined: Jul 2018
RE: (48G) Bernoulli numbers
Here is lua function for zigzag numbers, build 1 zigzag at a time

Code:
function zigzag(t)      -- Euler zigzag numbers
t = t or {[0] = 1}  -- t[0], t[-0] = T(0), T(1)
local i = #t + 1
t[i] = 0; t[-i] = 0 -- make room
for j=2-i, i    do t[j]=t[j]+t[j-1] end -- zig
for j=i-1,-i,-1 do t[j]=t[j]+t[j+1] end -- zag
return t            -- t[i], t[-i] = T(next even), T(next odd)
end

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

-94598037819122125295227433069493721872702841533066936133385696204311395415197247​711 / 33330
09-07-2023, 04:18 PM
Post: #6
 John Keith Senior Member Posts: 1,040 Joined: Dec 2013
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:

\<<
IF DUP 2 <                           @ If m < 2
THEN 1 SWAP
IF 0 SAME
THEN 1                             @ Return  1, 1 if m = 0
ELSE NEG 2                         @ Return -1, 2 if m = 1
END
ELSE
IF DUP 2 MOD                       @ If m odd,
THEN DROP 0 1                      @ return 0, 1
ELSE \-> m
\<< 1 2 m ^ 1 - 3 m 1 +          @ d = 1; n = 2^m-1
FOR k
IF DUP k MOD NOT             @ If k | 2^m-1
THEN k /
WHILE DUP k MOD NOT        @ Remove all factors of k from n
REPEAT k /
END
IF m k 1 - MOD NOT         @ If k-1 | m
THEN SWAP k * SWAP         @ d = d*k
END
END 2
STEP DROP
m ! 6.28318530718 m ^ / 2 *    @ 2 * m! / (2*pi)^m
m 2 / 2 MOD
IF NOT
THEN NEG                       @ Correct sign
END OVER *                     @ Multiply by d
m NEG 3 OVER ^ 2 ROT ^ + 1 + * @ Apply correction factor
FLOOR 2 * 1 +                  @ Numerator
SWAP 2 *                       @ Denominator = d * 2
\>>
END
END
\>>
09-08-2023, 08:09 PM
Post: #7
 John Keith Senior Member Posts: 1,040 Joined: Dec 2013
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
 Albert Chan Senior Member Posts: 2,705 Joined: Jul 2018
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,  711, 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
local n, d = 2^m-1, 1
for p= 3, sqrt(m)+1, 2 do           -- half p's
if n%p ~= 0 then continue end
repeat n=n/p until n%p ~= 0
if m%(p-1)==0 then d = d*p end  -- m = (p-1)*(q-1)
end
for p = 2, m do
local q = m/(p-1) + 1           -- half q's
if p >= q then break end        -- quit if p >= sqrt(m)+1
if m%(p-1)==0 and n%q==0 then d = d*q end
end
return 2*d
end

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
 Albert Chan Senior Member Posts: 2,705 Joined: Jul 2018
RE: (28 48) Bernoulli numbers
(09-09-2023 06:33 PM)Albert Chan Wrote:  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.

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
local n, d, hi = 2^m-1, 1, sqrt(m)
for p = 3, hi+1, 2 do               -- half p's
if n%p ~= 0 then continue end
repeat n = n/p until n%p ~= 0
if m%(p-1)==0 then d = d*p end  -- m = (p-1)*(q-1)
end
for k = 1, hi do
local q = m/k + 1               -- half q's
if m%k==0 and n%q==0 then d = d*q end
end
return 2*d
end

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
 John Keith Senior Member Posts: 1,040 Joined: Dec 2013
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
 Albert Chan Senior Member Posts: 2,705 Joined: Jul 2018
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
local n, d = 2^m-1, 1
for p = 3, m/2+1, 2 do
if n%p ~= 0 then continue end
repeat n = n/p until n%p ~= 0
if m%(p-1)==0 then d = d*p end
end
if n%(m+1)==0 then d = d*(m+1) end -- q for p=2
return 2*d
end
09-10-2023, 07:45 PM (This post was last modified: 09-10-2023 07:49 PM by John Keith.)
Post: #12
 John Keith Senior Member Posts: 1,040 Joined: Dec 2013
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:

\<< I\->R DUP 3.
IF <
THEN 1 -2 INV 6 INV 3. \->LIST SWAP 1. + GET
ELSE
IF DUP 2. MOD
THEN DROP 0
ELSE 2. / R\->I \-> n
\<< { 0 1 } 2 n
FOR k 0 :: + Scanr 0 :: + Scanl
NEXT Last NEWOB
n DUP I\->R 2. MOD :: NEG IFT *
4 n ^ DUP 1 - SWAP 2 / NEG * /
\>>
END
END
\>>

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:

\<< I\->R DUP 2. MOD - 4. MAX 2. / R\->I \-> n
\<< 1 -2 INV 6 INV { 0 1 } 2 n
FOR k 0 SWAP 0 :: + Scanr 0 :: + Scanl
DUP Last NEWOB
k I\->R 2. MOD :: NEG IFT
k * 4 k ^ DUP 1 - SWAP -2 / * / SWAP
NEXT DROP n 2 * 1 + \->LIST
\>>
\>>

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
 Gerald H Senior Member Posts: 1,627 Joined: May 2014
RE: (28 48 49 50) Bernoulli Numbers
On a 50g this revised programme

Code:
Size: 251.

CkSum: # 991Bh

::
CK1&Dispatch
# FF
::
ZINT 0
SWAP
FPTR2 ^DupQIsZero?
case2drop
ZINT 1
FPTR2 ^DupZIsOne?
case
::
ZINT 2
FPTR2 ^NDXFext
FPTR2 ^QSub
;
FPTR2 ^DupZIsEven?
NOT
case2drop
ZINT 0
DUP
ZINT 2
3UNROLL
FPTR2 ^Z2BIN
ZINT 0
ZINT 0
ZINT -1
4PICK
#2+
ONE_DO
ZINT 1
7PICK
INDEX@
FPTR2 ^#>Z
FPTR2 ^QDiv
FPTR2 ^QSub
FPTR2 ^QMul
5PICK
5PICK
FPTR2 ^PPow#
ROT
SWAP
5ROLL
ZINT 1
5UNROLL
2DUP
FPTR2 ^QMul
INDEX@
FPTR2 ^#>Z
FPTR2 ^QDiv
4ROLL
3UNROLL
LOOP
2DROP
4UNROLL3DROP
;
;

for input

102

returns

'3204019410860907078243020782116241775491817197152717450679002501086861530836678​158791/4326'

in

25.1814s
 « Next Oldest | Next Newest »

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