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

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。

这既不是 "overflow" 也不是 "underflow"，因为表示像 1e-308 这样小的数字并不困难。就是浮点运算时丢失有效位的问题。当你做1 + (1 / x)的时候，x越大，你加1的时候，(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|

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

避免大号加小号。在浮点加法a + b = a * (1 + b / a)中，如果b / a < 2.22e-16，那么我们a + b = a。这意味着在将多个正数相加时，从小到大依次累加更稳定。
避免用一个数减去相同大小的另一个数，否则您可能会得到 cancellation error。网页上有使用二次公式的经典例子

还建议您阅读 , a question asked a few days after your question. Using a series to approximate an irrational number is numerically stable as you won't get the absurd behavior seen in your question. But the finite number of valid significant digits imposes a different problem: numerical convergence, that is, you can only approximate the target value up to a certain number of significant digits. ，使用泰勒级数会在 19 项后收敛，因为 1 / factorial(19) 添加到 1 时在数值上已经是 0。

浮点数之间的乘法/除法不会导致有效数字的问题；它们可能会导致 "overflow" 或 "underflow"。然而，鉴于可表示的浮点值范围很广 (1e-308 ~ 1e+307)，"overflow" 和 "underflow" 应该很少见。真正的困难在于加法/减法，其中很容易丢失有效数字。参见 for an example on matrix computations. It is not impossible to get higher precision, but the work is probably more involved. For example, OP of the matrix example eventually used the GMP (GNU Multiple Precision Arithmetic Library) and associated R packages to proceed:

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"

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

R: approximating `e = exp(1)` using `(1 + 1 / n) ^ n` gives absurd result when `n` is large

r

exponential