(28 48 49 50) Bernoulli Numbers

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