如何修复自然常数e计算过程中的精度问题

原始代码 #

a = 1.0:0.1:[y/10 | (x,y) <- zip a (tail a)]
b = map (\n -> (1+n)**(1/n)) a

问题现象 #

前十几个计算结果能正常逼近e，但当取索引较大的元素时，结果出现异常：

Prelude> b!!2
2.7048138294215285
Prelude> b!!3
2.7169239322355936
Prelude> b!!4
2.7181459268249255
Prelude> b!!5
2.718268237192297
Prelude> b!!6
2.718280469095753
Prelude> b!!500
1.0

问题原因 #

当n趋近于极小值时，1+n会被Double浮点数近似为1.0——这是因为Double的精度有限（机器epsilon约为2.2e-16），无法区分1和1+极小值。此时(1.0)**(1/n)的结果恒为1.0，完全偏离了e的逼近值。

修复方案 #

改用数值稳定性更好的等价数学表达式：用指数和对数函数替代直接的幂运算。因为(1+n)^(1/n)等价于exp( log(1+n)/n )，这个形式在n极小时能更准确地计算逼近值，不会因1+n被舍入为1而失效。

修改后的代码：

a = 1.0:0.1:[y/10 | (x,y) <- zip a (tail a)]
b = map (\n -> exp (log(1+n)/n)) a

效果验证 #

修改后，即使取索引很大的元素，结果也能正常逼近e。例如b!!500会返回接近2.71828的数值，而非1.0。

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