Computes the Bernoulli numbers in the desired (binary) precision.
The computation happens via the `zeta`

function and the
formula

*B_k = -k ζ(1 - k),*

and hence the only non-zero odd Bernoulli number is *B_1 = +1/2*.
(Another tradition defines it, equally sensibly, as *-1/2*.)

1 | ```
Bernoulli(k, precBits = 128)
``` |

`k` |
non-negative integer vector |

`precBits` |
the precision in |

an `mpfr`

class vector of the same length as
`k`

, with i-th component the `k[i]`

-th Bernoulli number.

Martin Maechler

http://en.wikipedia.org/wiki/Bernoulli_number

`zeta`

is used to compute them.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 | ```
Bernoulli(0:10)
plot(as.numeric(Bernoulli(0:15)), type = "h")
curve(-x*zeta(1-x), -.2, 15.03, n=300,
main = expression(-x %.% zeta(1-x)))
legend("top", paste(c("even","odd "), "Bernoulli numbers"),
pch=c(1,3), col=2, pt.cex=2, inset=1/64)
abline(h=0,v=0, lty=3, col="gray")
k <- 0:15; k[1] <- 1e-4
points(k, -k*zeta(1-k), col=2, cex=2, pch=1+2*(k%%2))
## They pretty much explode for larger k :
k2 <- 2*(1:120)
plot(k2, abs(as.numeric(Bernoulli(k2))), log = "y")
title("Bernoulli numbers exponential growth")
Bernoulli(10000)# - 9.0494239636 * 10^27677
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.