R：如何计算 n-choose-k 中的大数？

Question

对于 class 作业，我需要创建一个函数来计算 n 选择 k。我就是这样做的，它适用于小数字（例如 6 选择 2），但我应该让它与 200 选择 50 一起工作，而它自然不会。答案太大，R 输出 NaN 或 Inf，表示：

> q5(200, 50)
[1] "NaN"
Warning message:
In factorial(n) : value out of range in 'gammafn'

我尝试使用对数和指数，但没有成功。

q5 <- function (n, k) {
  answer <- log(exp( factorial(n) / ( (factorial(k)) * (factorial(n - k)) )))
  paste0(answer)
}

Answer 1

大号套餐：

Brobdingnag“R 中非常大的数字”的包：
https://cran.r-project.org/web/packages/Brobdingnag/index.html
论文：https://www.researchgate.net/publication/251996764_Very_large_numbers_in_R_Introducing_package_Brobdingnag

library(Brobdingnag)
googol <- as.brob(10)^100 # googol:=10^100
googol
# [1] +exp(230.26) # exponential notation is convenient for very large numbers

gmp 多个精度算术包（大整数 和有理数、素数测试、矩阵计算）：
https://cran.r-project.org/web/packages/gmp/index.html

Answer 2

你也可以试试这个：

q5 <- function (n, k) {
  # nchoosek = (n-k+1)(n-k+2)...n / (1.2...k)
  return(prod(sapply(1:k, function(i)(n-k+i)/(i))))
}

q5(200, 50)
#[1] 4.538584e+47

或在日志域中

q5 <- function (n, k) {
  # ln (nchoosek) = ln(n-k+1) + ln(n-k+2) + ...+ ln(n) - ln(1) -ln(2) - ...- ln(k)
  return(exp(sum(sapply(1:k, function(i)(log(n-k+i) - log(i))))))
}
q5(200, 50)
#[1] 4.538584e+47

Answer 3

实际问题的答案是 R 无法显示它无法表示的数字，并且您的方程式中的某些项太大而无法表示。所以它失败了。然而，可以使用阶乘的近似值——它们使用对数，对数变大的速度要慢得多。

最著名的 Sterling 近似值不够准确，但 Ramanujan's approximation 派上了用场:)

ramanujan <- function(n){
  n*log(n) - n + log(n*(1 + 4*n*(1+2*n)))/6 + log(pi)/2
}

nchoosek <- function(n,k){
  factorial(n)/(factorial(k)*factorial(n-k))
} 

bignchoosek <- function(n,k){
  exp(ramanujan(n) - ramanujan(k) - ramanujan(n-k))
}

nchoosek(20,5)
# [1] 15504

bignchoosek(20,5)
# [1] 15504.06


bignchoosek(200,50)
# [1] 4.538584e+47

Answer 4

此解决方案计算帕斯卡三角形的完整行：

x <- 1
print(x)
for (i in 1:200) { x <- c(0, x) + c(x, 0); print(x) }
x[51]  ### 200 choose 50
## > x[51]
## [1] 4.538584e+47

（正如我为 How would you program Pascal's triangle in R? 提议的那样）
如果你想加速代码，那么不要 print(x)（输出是一个相对较慢的操作）。

把代码放在我们可以做的函数中

nchoosek <- function(n,k) {
    x <- 1
    for (i in 1:n) x <- c(0, x) + c(x, 0)
    x[k+1]  ### n choose k
} 

nchoosek(200, 50)    ### testing the function
## [1] 4.538584e+47

这是我的函数的更完善的版本：

nchoosek <- function(n, k) {
  if (k==0) return(1)
  if (k+k > n) k <- n-k
  if (k==0) return(1)
  x <- 1
  for (i in 1:k)     x <- c(0, x) + c(x, 0)
  for (i in 1:(n-k)) x <- x + c(0, head(x, -1))
  tail(x, 1)
} 
nchoosek(200, 50)    ### testing the function
## [1] 4.538584e+47

R：如何计算 n-choose-k 中的大数？

R: How can I calculate large numbers in n-choose-k?

r

factorial

exponent