R语言精度问题：为何1/0.00001运算会引发异常而其他量级不会？

现象复现 #

运行模拟时碰到一个奇怪的问题：直接传整数100000给rnorm()，返回的随机数长度是100000；但传1/0.00001时，返回长度变成了99999：

> length(rnorm(100000))
[1] 100000
> length(rnorm(1/0.00001))
[1] 99999

但换其他量级的运算就没问题：

> length(rnorm(1000000))
[1] 1000000
> length(rnorm(1/0.000001))
[1] 1000000

用Rmpfr包查看高精度计算结果，发现1/0.00001的实际存储值根本不是精确的100000：

> mpfr(1/0.0001,precBits = 53)
1 'mpfr' number of precision  53   bits 
[1] 10000
> mpfr(1/0.00001,precBits = 53)
1 'mpfr' number of precision  53   bits 
[1] 99999.999999999985
> mpfr(1/0.000001,precBits = 53)
1 'mpfr' number of precision  53   bits 
[1] 1000000
> mpfr(1/0.0000001,precBits = 53)
1 'mpfr' number of precision  53   bits 
[1] 10000000

问题根源 #

这完全是二进制浮点数的天生精度缺陷导致的：

• 计算机用二进制存储浮点数，但十进制的0.00001（也就是1e-5）没法被二进制浮点数精确表示，它的二进制近似值转成十进制后，会比真实的0.00001稍微小一点。
• 当你算1 / 0.00001时，实际上是用1除以这个略小的近似值，结果自然就会略小于100000（也就是上面看到的99999.999999999985）。
• rnorm()会自动把输入的数值转成整数，而且是直接截断小数部分（和as.integer()的逻辑一样），所以这个接近100000但略小的数，就被砍成了99999，生成的随机数长度也就跟着变了。

至于其他量级的数值（比如0.0001、0.000001），它们的二进制近似值在和1做除法后，结果的误差小到被浮点数精度当成了整数，所以转成整数时不会出问题。

解决办法 #

要避免这种坑，要么直接传整数给rnorm()，要么对计算结果做显式的四舍五入：

# 直接用整数，最稳妥
length(rnorm(100000))
# 显式四舍五入，修正精度误差
length(rnorm(round(1/0.00001)))

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