# computing Bernoulli numbers

Is there a good way to compute the ratio ( B[n] / n! ) that occurs so often in power series coefficients? Good in the sense that you get an answer that does not overflow a double; the largest n such that B[n] fits in a double is relatively small, but the ratio in question should never overflow, although I guess it will eventually underflow.

The Wikipedia article on Bernoulli numbers has a simple algorithm for the Bernoulli numbers, more than adequate for this range, but it's not obvious how to modify it to answer this question. To clarify, I'm not asking for an exact computation, but a decimal computation.

these are socalled "scaled" Bernoulli numbers; an efficient algorithm is discussed by Brent & Harvey, arXiv:1108.0286, Eq. 8.

alternatively, you can use any method that computes the Riemann zeta function (say via the Euler product), because of the identity

$$\frac{B_{2n}}{(2n)!}=(-1)^{n+1}2\zeta(2n)(2\pi)^{-2n}$$

(note that $B_{2n+1}=0$ for $n\geq 1$)

• I believe that the Brent/Harvey method is the one used in Sage for exact computation. I wasn't aware that they also gave an algorithm for the scaled Bernoulli numbers, so thanks for pointing that out. All the exact algorithms I found do not work well if adapted to fixed-precision computation (i.e. with doubles) because round-off error piles up and gives you wrong answers by n=54 or sooner. Commented May 25, 2013 at 14:04

Of course, We have several explicit and closed formulas which you can use for $\frac{B_n}{n!}$, for instance,

$\frac{B_m}{m} = \frac{1}{m}\sum_{k=1}^{m+1}\sum_{v=1}^{k}(-1)^{v+1}\binom{k-1}{v-1}\frac{v^m}k$

or

$\frac{B_k}{k}=\frac{1}{k}\sum_{n=0}^k \frac{1}{n+1}\sum_{j=0}^n(-1)^j \binom{n}{j}j^k$ ; $k>1$

So, you can get the desired result.Moreover you can study following paper for computing Bernoulli numbers , see