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

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

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

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

此解决方案将计算 roughlyPI 1000，但由于 2 ^ 1024 :: Double == Infinity，它在 1023 上以 NaN 失败。注意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)

我们可以做得比在 1000 后爆破 Double 更好，方法是使用 Rationals 进行计算，然后使用 realToFrac 转换为 Double（归功于 @leftaroundabout）：

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)

更多参考请见Wikipedia page on approximations of PI

附：抱歉，公式太长了，stackoverflow 不支持 LaTex

Answer 2

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

*Main> roughlyPI 91
3.1415926535897922

如前所述，问题在于，当您尝试使近似值更好时，阶乘项变得太大而无法用双精度浮点数表示。解决这个问题的最简单（尽管有点蛮力）的方法是用 rational 算术代替所有计算。因为 Haskell 中的数值运算是多态的，所以它可以使用与您几乎相同的代码，只有 ** 运算符不能使用，因为它允许小数指数（通常是无理数）。相反，您应该使用 integer 指数，这在概念上是正确的。这需要几个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

Answer 3

在这种情况下要走的路是计算连续项的比率，并通过滚动比率的乘法来计算项：

-- 1. -------------
pi1 n = Sum { k = 0 .. n } T(k)
  where
  T(k) = 2^(k+1)(k!)^2 / (2k+1)!

-- 2. -------------
ts2 = [ 2^(k+1)*(k!)^2 / (2k+1)! | k <- [0..] ]
pis2 = scanl1 (+) ts2
pi2 n = pis2 !! n

-- 3. -------------
T(k)  = 2^(k+1)(k!)^2 / (2k+1)!
T(k+1) = 2^(k+2)((k+1)!)^2 / (2(k+1)+1)!
       = T(k) 2 (k+1)^2 / (2k+2) (2k+3)
       = T(k)   (k+1)^2 / ( k+1) (2k+3)
       = T(k)   (k+1)   /        (k+1 + k+2)
       = T(k)           /        (1 + (k+2)/(k+1))
       = T(k)           /        (2 + 1    /(k+1))

-- 4. -------------
ts4 = scanl (/) 2 [ 2 + 1/(k+1) | k <- [0..]] :: [Double]
pis4 = scanl1 (+) ts4
pi4 n = pis4 !! n

通过这种方式，我们尽可能地共享和重用计算。这导致最有效的代码，希望导致最小的累积数值错误。这个公式也非常简单，甚至可以进一步简化为ts5 = scanl (/) 2 [ 2 + recip k | k <- [1..]]。

试一试：

> pis2 = scanl1 (+) $ [ fromIntegral (2^(k+1))*fromIntegral (product[1..k])^2 / 
                         fromIntegral (product[1..(2*k+1)]) | k <- [0..] ] :: [Double]

> take 8 $ drop 30 pis2
[3.1415926533011587,3.141592653447635,3.141592653519746,3.1415926535552634,
 3.141592653572765,3.1415926535813923,3.141592653585647,3.141592653587746]

> take 8 $ drop 90 pis2
[3.1415926535897922,3.1415926535897922,NaN,NaN,NaN,NaN,NaN,NaN]

> take 8 $ drop 30 pis4
[3.1415926533011587,3.141592653447635,3.141592653519746,3.1415926535552634,
 3.141592653572765,3.1415926535813923,3.141592653585647,3.141592653587746]

> take 8 $ drop 90 pis4
[3.1415926535897922,3.1415926535897922,3.1415926535897922,3.1415926535897922,
 3.1415926535897922,3.1415926535897922,3.1415926535897922,3.1415926535897922]

> pis4 !! 1000
3.1415926535897922