R: 当`n`非常大时，使用`(1 + 1 / n) ^ n`来近似计算`e = exp(1)`会得到荒谬的结果。

Question

所以，我刚刚在R中手动计算e的值，发现了一些让我有点不安的事情。

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

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

嗯，这不对。发生了什么？我是否溢出了数据类型，需要使用特定的包？如果是这样，当您溢出数据类型时没有警告吗？

- self_learnt

---

怎么样？e是我输入的函数的水平渐近线。它不应该接近1，而应该接近e。 - self_learnt

@self_learnt 也许你可以看看这个问题：为什么这些数字不相等？ - Ronak Shah

Answer 1

你的问题在于机器精度。一旦 (1 / x) < 2.22e-16，1 + (1 / x) 就只是 1。在有限精度数值计算中，数学极限会崩溃。你问题中的最终 x 已经是 5e+15，非常接近这个边界。尝试 x <- x * 10，你的 y 将会是 1。

这既不是 "overflow" 也不是 "underflow"，因为表示一个小到 1e-308 的数字并没有困难。这是浮点运算期间有效数字丢失的问题。当你执行 1 + (1 / x) 时，x 越大，(1 / x) 部分可以保留的有效数字就越少，最终你会完全失去那个 (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|

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

- 避免将一个大数和一个小数相加。在浮点数加法中，如果 b / a < 2.22e-16，那么 a + b = a，即当累加一系列正数时，从小到大逐个相加更加稳定。 - 避免从同等大小的数字中减去另一个数字，否则可能会出现cancellation error。该网页上有一个使用二次公式的经典例子。

建议您阅读在50次迭代后，常数“π”的近似值不会更好, 这是在您提问几天后发布的一个问题。使用级数来逼近无理数是数值上稳定的，因为您不会遇到在您的问题中看到的荒谬行为。但是，有效数字的有限数量带来了不同的问题：数值收敛，也就是说，您只能将目标值近似到一定数量的有效数字。MichaelChirico的答案使用泰勒级数，在19项之后收敛，因为当加到1时，1 / factorial(19)已经是数值上的0了。

浮点数之间的乘除不会影响有效数字，但可能会导致"溢出"或"下溢"。然而，由于可表示的浮点数值范围很大（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?

- Zheyuan Li

Answer 2

你也可以尝试使用泰勒级数逼近 exp(1)，即

e^x = \sum_{k = 0}{\infty} x^k / k!

因此，我们可以通过截断这个和来近似得到e = e^1；在 R 中：

sprintf('%.20f', exp(1))
# [1] "2.71828182845904509080"
sprintf('%.20f', sum(1/factorial(0:10)))
# [1] "2.71828180114638451315"
sprintf('%.20f', sum(1/factorial(0:100)))
# [1] "2.71828182845904509080"

- MichaelChirico

---

@李哲源 是的，实际上 sprintf('%.70f', sum(1/factorial(0:17))) 从 R 的角度来看已经数值等同于 exp(1)。 - MichaelChirico

所以，最好写成sum(1/factorial(17:10))？有趣！ - MichaelChirico