library(Rmpfr)

mp <- function(x) mpfr(x, precBits=100)

rel.err <- function(x, y)
{
	if (any(sign(x) != sign(y)))
		return(Inf)
	ind <- x != 0 | y != 0
	x <- abs(x[ind])
	y <- abs(y[ind])
	max(abs(x - y)/pmin(x, y), 0)
}

m <- 1000
n <- round(runif(m, min = -1000, max = 1200)) + (2*runif(m) - 1)^7 # concentrates near integers
mpfr.n <- mp(n)
mpfr.x <- mp(rep(1, times=length(n))) # initialize Rmpfr computed choose(n, k)
max.err <- mp(0) # initialize maximum relative error
for (k in 0:209)
{
	y <- choose(n, k)
	max.err <- max(max.err, rel.err(mpfr.x, y))
	if (k %% 10 == 9)
		cat("k <=", k, " max rel err =", as.double(max.err), "\n")
	mpfr.x <- mpfr.x*(mpfr.n - k)/(k + 1) # prepare for the next k
}