(28 48 49 50) Bernoulli Numbers

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