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