R中用星与杠逻辑实现骰子和分布时15枚以上骰子结果错误问题

• 第一：代码硬编码了骰子面数为6，与函数的sides参数不匹配：公式中写死了6*j，如果后续修改sides参数会直接出错，即使默认sides=6也属于不规范写法，不过这不是骰子≥15出错的直接原因。
• 第二：判断条件不完整：当前代码仅判断n > -1就计算组合数，虽然R的choose函数在k<0或n<k时会返回0，但显式增加k >= 0和n >= k的判断可以避免不必要的计算，也能规避部分边界场景的异常。
• 第三：核心错误是浮点数精度溢出：R默认使用双精度浮点数存储数值，双精度浮点数仅能精确表示≤2^{53（约9e15）的整数。当骰子数量≥15时，部分和对应的组合数已经接近/超过32位整数的最大值（2}31-1≈2e9），如果使用32位R会直接出现精度丢失；当骰子数量≥21时，总可能性6^{21≈2.19e16超过2}53，即使64位R也会出现双精度浮点数无法精确表示整数的问题，导致结果错误。

1. 替换默认的整数/浮点数运算为高精度整数运算，使用R的gmp包提供的大整数运算函数，避免精度丢失。
2. 把硬编码的6替换为sides参数，保证函数通用性。
3. 完善边界判断条件，减少无效计算。

修复后代码 #

# 需先安装gmp包：install.packages("gmp")
library(gmp)

dice_sum_distribution_sb <- function(dice = 2, sides = 6){
  # 总可能性用大整数计算
  all <- as.bigz(sides)^dice
  
  minSum <- dice
  maxSum <- dice * sides
  sumRange <- minSum:maxSum
  sumCounts <- rep(as.bigz(0), length(sumRange))
  
  for(i in seq_along(sumRange)){
    value <- sumRange[i]
    count <- as.bigz(0)
    # j的上限直接取到最大有效值，避免不必要的循环
    max_j <- (value - dice) %/% sides
    for(j in 0:max_j){
      n <- value - 1 - sides * j
      k <- value - dice - sides * j
      # 仅当组合数有效时计算
      if(n >= k && k >= 0){
        term <- (-1)^j * chooseZ(dice, j) * chooseZ(n, k)
        count <- count + term
      }
    }
    sumCounts[i] <- count
  }
  
  distribution <- data.frame(
    sum = sumRange,
    sumCount = as.numeric(sumCounts), # 数值超过2^53时不要转numeric，直接保留bigz类型
    pdf = as.numeric(sumCounts / all)
  )
  return(distribution)
}

骰子数量≤20时，修复后的代码输出结果与暴力枚举/动态规划结果完全一致；骰子数量>20时，如果需要精确的组合数，不要把sumCount转为numeric，直接使用bigz类型存储即可。

内容的提问来源于stack exchange，提问作者wrudi