On Rosetta Code, the list of Tasks uncompleted in R contains the following challenge: print the Bernoulli numbers from B0 to B60. After seeing some solutions, including one in in Python , I wanted to try in R. There are examples of code on that website, and a nice description of the algorithm here.
However, despite many different approaches, I can’t seem to get any answers past 17, apparently because the numbers involved become too large! Any thoughts?
this code available as a gist
In addition to the matrix approach below, I tried using lists
"%gcd%" <- function(u, v) {ifelse(u %% v != 0, v %gcd% (u%%v), v)}
n <- 20L
js <- 0:n
Jp1 <- js + 1L
Ns <- matrix(NA,nrow = n + 1, ncol = n + 1)
Ds <- Ns
Ns[1,] <- 1L
Ds[1,] <- js+1L
for(i in 2:(n + 1)){
Ns[i,] <- Ns[i-1,]*Ds[i-1,c(2:(n + 1),NA)] - Ds[i-1,]*Ns[i-1,c(2:(n + 1),NA)]
Ns[i,] <- Ns[i,]*Jp1
Ds[i,] <- Ds[i-1,]*Ds[i-1,c(2:(n + 1),NA)]
gcds <- Ns[i,] %gcd% Ds[i,]
Ns[i,] <- Ns[i,] %/% gcds
Ds[i,] <- Ds[i,] %/% gcds
}## Warning: NAs produced by integer overflow
## Warning: NAs produced by integer overflow
## Warning: NAs produced by integer overflow
## Warning: NAs produced by integer overflow
## Warning: NAs produced by integer overflowresult <- data.frame(n = js,Bn = paste0(Ns[,1],"/",Ds[,1]))
library(xtable)
print(xtable(result),include.rownames = FALSE, type = "html")| n | Bn |
|---|---|
| 0 | 1/1 |
| 1 | 1/2 |
| 2 | 1/6 |
| 3 | 0/1 |
| 4 | -1/30 |
| 5 | 0/1 |
| 6 | 1/42 |
| 7 | 0/1 |
| 8 | -1/30 |
| 9 | 0/1 |
| 10 | 5/66 |
| 11 | 0/1 |
| 12 | -691/2730 |
| 13 | 0/1 |
| 14 | 7/6 |
| 15 | 0/1 |
| 16 | -3617/510 |
| 17 | NA/NA |
| 18 | NA/NA |
| 19 | NA/NA |
| 20 | NA/NA |
So then I tried using Rmpfr, which allows for more precise integers. It works much better, but still doesn’t get as far as 60:
library(Rmpfr)Loading C code of R package ‘Rmpfr’: GMP using 64 bits per limb
library(MASS)
n <- 40L
js <- 0:n
Jp1 <- js + 1L
Ns <- matrix(NA,nrow = n + 1, ncol = n + 1)
Ds <- Ns
Ns[1,] <- 1L
Ds[1,] <- js+1L
Ns <- mpfr(Ns,precBits = 200)
Ds <- mpfr(Ds,precBits = 200)
for(i in 2:(n + 1)){
Ns[i,-(n + 1)] <- Ns[i-1,-(n + 1)]*Ds[i-1,2:(n + 1)] - Ds[i-1,-(n + 1)]*Ns[i-1,2:(n + 1)]
Ns[i,] <- Ns[i,]*Jp1
Ds[i,-(n + 1)] <- Ds[i-1,-(n + 1)]*Ds[i-1,2:(n + 1)]
}
result <- mpfrArray(0,precBits = 200, dim = c(nrow(Ns),3))
result[,1] <- Ns[,1]
result[,2] <- Ds[,1]
result[,3] <- result[,1] / result[,2]
nums <- asNumeric(result[,3])
fracs <- fractions(nums)
result2 <- data.frame(B = js,Bn = attr(fracs,"frac"))
print(xtable(result2),include.rownames = FALSE, type = "html")| B | Bn |
|---|---|
| 0 | 1 |
| 1 | 1/2 |
| 2 | 1/6 |
| 3 | 0 |
| 4 | -1/30 |
| 5 | 0 |
| 6 | 1/42 |
| 7 | 0 |
| 8 | -1/30 |
| 9 | 0 |
| 10 | 5/66 |
| 11 | 0 |
| 12 | -691/2730 |
| 13 | 0 |
| 14 | 7/6 |
| 15 | 0 |
| 16 | -3617/510 |
| 17 | 0 |
| 18 | 43867/798 |
| 19 | 0 |
| 20 | -174611/330 |
| 21 | 0 |
| 22 | 854513/138 |
| 23 | 0 |
| 24 | -236364091/2730 |
| 25 | 0 |
| 26 | 8553103/6 |
| 27 | 0 |
| 28 | -58466341266499/2141763 |
| 29 | NaN |
| 30 | NaN |
| 31 | NaN |
| 32 | NaN |
| 33 | NaN |
| 34 | NaN |
| 35 | NaN |
| 36 | NaN |
| 37 | NaN |
| 38 | NaN |
| 39 | NaN |
| 40 | NaN |