对无限级数的有限前缀求和

Question

该数字?可以通过以下无穷级数和来计算：

我想定义一个 Haskell 函数roughlyPI，给定一个自然数k，计算从0到k值的序列总和。

示例：（roughlyPi 1000或其他）=>3.1415926535897922

我所做的是（在 VS Code 中）：

roughlyPI :: Double -> Double 
roughlyPI 0 = 2 
roughlyPI n = e1/e2 + (roughlyPI (n-1)) 
         where 
             e1 = 2**(n+1)*(factorial n)**2 
             e2 = factorial (2*n +1) 
             factorial 0 = 1 
             factorial n = n * factorial (n-1)

但它并没有真正起作用......

*Main> roughlyPI 100 
NaN

我不知道怎么了。顺便说一下，我是 Haskell 的新手。

我真正想要的是能够输入一个PI最后会给我的数字。不可能那么难...

Answer 1

正如评论中提到的，我们需要避免大的划分，而是在阶乘中散布较小的划分。我们Double用于表示 PI 但也Double有其局限性。例如1 / 0 == Infinity和(1 / 0) / (1 / 0) == Infinity / Infinity == NaN。

幸运的是，我们可以使用代数来简化公式，并希望延迟我们的Doubles的爆发。通过在我们的阶乘中进行除法，数字不会变得太笨重太快。

此解决方案将计算roughlyPI 1000，但在 1023 上失败，NaN因为2 ^ 1024 :: Double == Infinity。请注意 的每次迭代如何fac进行除法和乘法以帮助防止数字爆炸。如果你试图用计算机来近似 PI，我相信有更好的算法，但我试图让它在概念上尽可能接近你的尝试。

roughlyPI :: Integer -> Double
roughlyPI 0 = 2
roughlyPI k = e + roughlyPI (k - 1)
  where
    k' = fromIntegral k
    e = 2 ** (k' + 1) * fac k / (2 * k' + 1)
      where
        fac 1 = 1 / (k' + 1)
        fac p = (fromIntegral p / (k' + fromIntegral p)) * fac (p - 1)

Double通过进行计算Rationals然后转换为Doublewith realToFrac（归功于@leftaroundabout），我们可以做得比在 1000 之后爆炸更好：

roughlyPI' :: Integer -> Double
roughlyPI' = realToFrac . go
  where
    go 0 = 2
    go k = e + go (k - 1)
      where
        e = 2 ^ (k + 1) * fac k / (2 * fromIntegral k + 1)
          where
            fac 1 = 1 % (k + 1)
            fac p = (p % (k + p)) * fac (p - 1)

如需进一步参考，请参阅有关 PI 近似值的维基百科页面

PS 抱歉，公式庞大，stackoverflow 不支持 LaTex

Answer 2

首先请注意，您的代码实际上有效：

*Main> roughlyPI 91
3.1415926535897922

正如已经说过的，问题是当您尝试使近似更好时，阶乘项变得太大而无法用双精度浮点数表示。解决这个问题的最简单的——尽管有点暴力——的方法是用有理算术来代替所有的计算。因为 Haskell 中的数值运算是多态的，所以它可以使用与您几乎相同的代码，只有**运算符不能使用，因为它允许分数指数（通常是无理数）。相反，您应该使用整数指数，这在概念上是正确的。这需要一些fromIntegral：

roughlyPI :: Integer -> Rational 
roughlyPI 0 = 2 
roughlyPI n = e1/e2 + (roughlyPI (n-1)) 
         where 
             e1 = 2^(n+1)*fromIntegral (factorial n^2)
             e2 = fromIntegral . factorial $ 2*n + 1
             factorial 0 = 1 
             factorial n = n * factorial (n-1)

这现在也适用于更高程度的近似，尽管需要很长时间来处理所涉及的巨大分数：

*Main> realToFrac $ roughlyPI 1000
3.141592653589793