collect.errors <- function(err, x, y)
{
    min.xy <- min(abs(x), abs(y))
    if (min.xy == 0) {
        tab.err[1, 2] <- max(tab.err[1, 2], abs(x - y))
        return(tab.err)
    }
    else {
        for (i in 2:nrow(tab.err))
            if (min.xy <= tab.err[i, 1]) {
                err <- max(x/y, y/x)
                if (err < 0) err <- Inf
                tab.err[i, 2] <- max(tab.err[i, 2], err - 1)
                return(tab.err)
            }
    }
}

options(warn=2)

seed <- 2
set.seed(seed)

ni <- 0:1000
m <- 10000
nr <- floor(runif(m, min=-50, max=150)) + (2*runif(m) - 1)^7 # concentrates near integers
min.err <- Inf # minimum inexact result for integer n
# table of maximum error for real n in different ranges
tab.err <- rbind(c(bound=0,error=0), c(1e-20,0), c(1e-15,0), c(1e-10,0), c(1e-5,0), c(1,0), c(Inf, 0))
xi <- rep(1, times=length(ni)) # exact values (Pascal triangle)
xr <- rep(1, times=length(nr)) # R level computed values
cases.warning <- NULL
for (k in 0:50) {
    print(k)
    # check integer n
    yi <- choose(ni + k, k)
    min.err <- min(min.err, xi[xi != yi])
    xi <- cumsum(xi) # prepare for the next k
    # check real n
    for (i in 1:length(nr)) {
        yr <- try( choose(nr[i], k) )
        if (class(yr) == "numeric") {
            tab.err <- collect.errors(tab.err, xr[i], yr)
        } else {
            cases.warning <- rbind(cases.warning, c(i, k))
        }
    }
    xr <- xr*(nr - k)/(k + 1) # prepare for the next k
}

cat("Cases with warning\n")
print(cases.warning)
cat("minimum log2() of a wrong result for integer n :", log2(min.err), "\n")
cat("maximum error for real n :\n")
print(tab.err)