haskell 纳皮耶常数计算精度的确定

yptwkmov  于 2023-03-08  发布在  其他
关注(0)|答案(1)|浏览(162)

我试着用无限数组来计算纳皮耶常数 e

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

数字的前一打的结果是可以的,但不是更小的数字;

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

你知道怎么修吗?

qv7cva1a

qv7cva1a1#

你的方法相当于更容易理解的

b = [ (1+1/n)^n | i<-[0..], let n=10^i ]

请注意,每一项都是单独计算的,因此将它们作为一个无限列表来计算没有多大意义--你还不如直接计算具有足够准确度的一项,而跳过所有其他项。
你观察到的问题是这个算法在数值上是不稳定的;它不应该在实际中使用。正如chi所指出的,问题在于,在浮点运算中,1+n最终的计算结果恰好是1,并且n足够小。1的任意次幂总是1。
你可能会认为这可以通过切换到精确算术而不是浮点来规避,实际上,理论上你可以用Rational类型来计算上面的值,并且会正确地收敛,但是你会发现这是完全不可行的,因为当分数被提升到一个幂的时候,分数会变得越来越大,并且在你获得任何优于浮点版本的优势之前,计算就已经停止了。
此外,这只是因为我将它从n=10^^(-i)**(1/n)改为整数n=10^i,这允许使用^运算符。**支持小数参数,但对于有理数来说这是不可能的(小数指数对应于根,通常是无理数)。
真正具有讽刺意味的是,**是按照 * 指数 * 实现的:

x**y = exp(log x*y)

所以你的原始版本实际上计算为

exp(log(1 + n)/n)
  • 同样,不稳定(你得到log(1+n) == 0),但是在这一点上,仅仅用解析的方法取极限是微不足道的,也就是
exp((0 + n + 𝓞(n²))/n)
 → exp(n/n)
 ≡ exp 1
  • 当然这是欧拉/纳皮耶常数...但它是自我参照的,不是吗?
    如果你想以一种真正有意义的方式计算这个常数,你应该使用级数展开。

相关问题