Bernoulli_Matrix.R

## via base ####

"%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
}

result <- data.frame(n = js,Bn = paste0(Ns[,1],"/",Ds[,1]))


## using Rmpfr ####


library(Rmpfr)
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"))