R:当 `e = exp(1)` 很大时，使用 `(1 + 1/n) ^ n` 近似 `n` 会产生荒谬的结果

Question

所以，我只是在玩手动计算 R 中 e 的值，我注意到一些让我有点不安的事情。

e 的值使用 R 的 exp() 命令...

exp(1)
#[1] 2.718282

现在，我将尝试使用 x = 10000

手动计算它

x <- 10000
y <- (1 + (1 / x)) ^ x
y
#[1] 2.718146

不完全是，但我们会尝试使用 x = 100000

x <- 100000
y <- (1 + (1 / x)) ^ x
y
#[1] 2.718268

比较暖和，但还是有点偏...

x <- 1000000
y <- (1 + (1 / x)) ^ x
y
#[1] 2.71828

现在，让我们用一个巨大的来试试吧

x <- 5000000000000000
y <- (1 + (1 / x)) ^ x
y
#[1] 3.035035

嗯，这是不对的。这里发生了什么？我是否溢出了数据类型并需要使用某个包来代替？如果是这样，当您溢出数据类型时是否没有警告？

Answer 1

您的机器精度有问题。尽快(1 / x) < 2.22e-16 , 1 + (1 / x)只是 1。数学极限在有限精度数值计算中被打破。您的最终 x问题中已经是5e+15 ，非常接近这个边缘。试试 x <- x * 10 ，以及您的 y将是 1 .

这既不是“溢出”也不是“下溢”，因为表示像 1e-308 这样小的数字并不困难。 .这是浮点运算过程中丢失有效数字的问题。当你这样做 1 + (1 / x) , 越大 x是，(1 / x) 中的有效数字较少。部分加到 1 时可以保留，最终丢失 (1 / x)完全术语。

               ## valid 16 significant digits
1 + 1.23e-01 = 1.123000000000000|
1 + 1.23e-02 = 1.012300000000000|
   ...            ...
1 + 1.23e-15 = 1.000000000000001|
1 + 1.23e-16 = 1.000000000000000|

任何数值分析书都会告诉你以下内容。

• 避免添加一个大数字和一个小数字。浮点加法 a + b = a * (1 + b / a) , 如果 b / a < 2.22e-16 , 那里有我们 a + b = a .这意味着当将多个正数相加时，将它们从最小到最大累加更稳定。
• 避免从一个相同大小的数字中减去一个数字，否则您可能会得到 cancellation error .该网页有一个使用二次公式的经典示例。

还建议您阅读 Approximation to constant "pi" does not get any better after 50 iterations ，在您提出问题几天后提出的问题。使用 series近似一个无理数在数值上是稳定的，因为您不会在问题中看到荒谬的行为。但是有效有效数字的有限数量带来了一个不同的问题:数值收敛，也就是说，您只能将目标值逼近到一定数量的有效数字。 MichaelChirico's answer使用泰勒级数将在 19 项后收敛，因为 1 / factorial(19)加到 1 时已经是数字 0。

float 之间的乘法/除法不会导致有效数字出现问题；它们可能会导致“溢出”或“下溢”。但是，鉴于可表示的浮点值范围很广(1e-308 ~ 1e+307)，“上溢”和“下溢”应该很少见。真正的困难在于加法/减法，其中有效数字很容易丢失。见 Can I stably invert a Vandermonde matrix with many small values in R?有关矩阵计算的示例。获得更高的精度并非不可能，但工作可能更复杂。例如，矩阵示例的 OP 最终使用了 GMP (GNU Multiple Precision Arithmetic Library)和相关的 R 包继续:How to put Rmpfr values into a function in R?