Summing a finite prefix of an infinite series

user_14786876 0
Viewed 187

The number π can be calculated with the following infinite series sum:

Formula

I want to define a Haskell function roughlyPI that, given a natural number k, calculates the series sum from 0 to the k value.

Example: roughlyPi 1000 (or whatever) => 3.1415926535897922

What I did was this (in 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)

but it doesn't really work....

*Main> roughlyPI 100 
NaN

I don't know what's wrong. I'm new to Haskell, by the way.

All I really want is to be able to type in a number that will give me PI at the end. It can't be that hard...

numerical-methods numeric precision functional-programming haskell
3 Answers

As mentioned in the comments, we need to avoid large divisions and instead intersperse smaller divisions within the factorials. We use Double for representing PI but even Double has its limits. For instance 1 / 0 == Infinity and (1 / 0) / (1 / 0) == Infinity / Infinity == NaN.

Luckily, we can use algebra to simplify the formula and hopefully delay the blowup of our Doubles. By dividing within our factorial the numbers don't grow too unwieldy too quickly.

f1

f2

f3

This solution will calculate roughlyPI 1000, but it fails on 1023 with NaN because 2 ^ 1024 :: Double == Infinity. Note how each iteration of fac has a division as well as a multiplication to help keep the numbers from blowing up. If you are trying to approximate PI with a computer, I believe there are better algorithms, but I tried to keep it as conceptually close to your attempt as possible.

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)

We can do better than having a blowup of Double after 1000 by doing computations with Rationals then converting to Double with realToFrac (credit to @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)

For further reference see Wikipedia page on approximations of PI

P.S. Sorry for the bulky equations, stackoverflow does not support LaTex

user_40510200

First note that your code actually works:

*Main> roughlyPI 91
3.1415926535897922

The problem, as was already said, is that when you try to make the approximation better, the factorial terms become too big to be representable in double-precision floats. The simplest – albeit somewhat brute-force – way to fix that is to do all the computation in rational arithmetic instead. Because numerical operations in Haskell are polymorphic, this works with almost the same code as you have, only the ** operator can't be used since that allows fractional exponents (which are in general irrational). Instead, you should use integer exponents, which is anyway the conceptually right thing. That requires a few 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)

This now works also for much higher degrees of approximation, although it takes a long time to carry around the giant fractions involved:

*Main> realToFrac $ roughlyPI 1000
3.141592653589793
user_7459030

The way to go in such cases is to calculate the ratio of consecutive terms and calculate the terms by rolling multiplications of the ratios:

-- 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

This way we share and reuse the calculations as much as possible. This leads to the most efficient code, hopefully leading to the smallest cumulative numerical error. The formula also turned out to be exceptionally simple, and could even be simplified further as ts5 = scanl (/) 2 [ 2 + recip k | k <- [1..]].

Trying it out:

> 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
user_8498910
Related