(28 48 49 50) Bernoulli Numbers

[John Keith]

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
  \>>
\>>

[Gerald H]

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
    FPTR2 ^RADDext
    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 ^QAdd
    FPTR2 ^CASCOMPEVAL
    SWAP
    5ROLL
    ZINT 1
    FPTR2 ^QAdd
    5UNROLL
    2DUP
    FPTR2 ^QMul
    INDEX@
    FPTR2 ^#>Z
    FPTR2 ^QDiv
    4ROLL
    FPTR2 ^QAdd
    FPTR2 ^CASCOMPEVAL
    3UNROLL
    LOOP
    2DROP
    4UNROLL3DROP
  ;
;

for input

1
100

returns

'
-94598037819122125295227433069493721872702841533066936133385696204311395415197247​711
/33330'

in

42.2994s

[John Keith]

Wow, that's fast! The built-in command takes about 77 seconds on my 50g. What algorithm are you using?

[Gerald H]

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

[Albert Chan]

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

[John Keith]

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
\>>

[John Keith]

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.

[Albert Chan]

(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
    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

[Albert Chan]

(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

[John Keith]

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.

[Albert Chan]

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

[John Keith]

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
  \>>
\>>

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.

[Gerald H]

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
    FPTR2 ^RADDext
    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
    FPTR2 ^QAdd
    SWAP
    5ROLL
    ZINT 1
    FPTR2 ^QAdd
    5UNROLL
    2DUP
    FPTR2 ^QMul
    INDEX@
    FPTR2 ^#>Z
    FPTR2 ^QDiv
    4ROLL
    FPTR2 ^QAdd
    3UNROLL
    LOOP
    2DROP
    4UNROLL3DROP
  ;
;

for input

102

returns

'3204019410860907078243020782116241775491817197152717450679002501086861530836678​158791/4326'

in

25.1814s