数字接近1的幂次方的强大威力

Question

我猜测有一些标准的技巧我没能找到：无论如何，我想以数值稳定的方式计算一个非常接近1（比如1-p，其中p < 1e-17）的大数次幂。在我的系统上，1-p被截断为1。

使用对数的泰勒展开式，我得到以下界限：

exp(-np-np^2/2) ≤ exp(n*log(1-p)) ≤ exp(-np)

是否有更聪明的方法？

- Tobias Madsen

---

一个替代方案是，由于当x变大时(1-x^-1)^x趋近于1/e，因此您可以将其计算为1/e的幂。 - James

@James 这不会导致与我的上限近似相同的结果吗？也就是说，(1-p)^n = [[1-1/(1/p)]^(1/p)]^(np) 约等于 (1/e)^np。 - Tobias Madsen

啊，是的。如果p很小，p的平方非常小，因此我认为您可以安全地使用exp（-np）。下限可以重写为exp（-np）^（1+p/2）。 - James

Answer 1

通过使用 log1p 函数，您可以更准确地计算 |x| <= 1 的 log(1+x)。

例如：

> p <- 1e-17
> log(1-p)
[1] 0
> log1p(-p)
[1] -1e-17

还有一个：

> print((1+1e-17)^100, digits=22)
[1] 1
> print(exp(100*log1p(-1e-17)), digits=22)
[1] 0.9999999999999990007993

然而，在这里，我们受到基于double类型的FP算术准确性的限制（请参见计算机科学家应知道的浮点算术）。

另一种方法是使用例如Rmpfr（也称为多精度浮点可靠）软件包：

> options(digits=22)
> library(Rmpfr)
> .N <- function(.) mpfr(., precBits = 200) # see the package's vignette
> (1-.N(1e-20))^100
1 'mpfr' number of precision  200   bits 
[1] 0.99999999999999999900000000000000005534172854579042829381053529

该软件包使用gsl和mpfr库来实现任意精度的浮点数运算（当然会导致计算速度变慢）。

- gagolews

---

很酷，我不知道log1p。如果有一些代数解决方法，不需要任意精度，那就太好了。 - Tobias Madsen

考虑到这篇论文，我认为很难找到任何解决方法... 这是一种普遍的计算限制。 - gagolews