Almkvist Berndt 1988 begins with an investigation of why the agm is such an efficient algorithm, and proves that it converges quadratically. This is an efficient method to calculate 



π


{\displaystyle \pi }

.

With the same notations used in Arithmetic-geometric mean, we can summarize the paper by writing:





π
=



4


a
g
m

(
1
,
1

/



2



)

2




1
−

∑

n
=
1


∞



2

n
+
1


(

a

n


2


−

g

n


2


)





{\displaystyle \pi ={\frac {4\;\mathrm {agm} (1,1/{\sqrt {2}})^{2}}{1-\sum \limits _{n=1}^{\infty }2^{n+1}(a_{n}^{2}-g_{n}^{2})}}}



This allows you to make the approximation, for any large   N:





π
≈



4


a

N


2




1
−

∑

k
=
1


N



2

k
+
1


(

a

k


2


−

g

k


2


)





{\displaystyle \pi \approx {\frac {4\;a_{N}^{2}}{1-\sum \limits _{k=1}^{N}2^{k+1}(a_{k}^{2}-g_{k}^{2})}}}



The purpose of this task is to demonstrate how to use this approximation in order to compute a large number of decimals of 



π


{\displaystyle \pi }

.

import Prelude hiding (pi)
import Data.Number.MPFR hiding (sqrt, pi, div)
import Data.Number.MPFR.Instances.Near ()
 
-- A generous overshoot of the number of bits needed for a
-- given number of digits.
digitBits :: (Integral a, Num a) => a -> a
digitBits n = (n + 1) `div` 2 * 8
 
-- Calculate pi accurate to a given number of digits.
pi :: Integer -> MPFR
pi digits =
  let eps = fromString ("1e-" ++ show digits)
            (fromInteger $ digitBits digits) 0
      two = fromInt Near (getPrec eps) 2
      twoi = 2 :: Int
      twoI = 2 :: Integer
      pis a g s n =
        let aB = (a + g) / two
            gB = sqrt (a * g)
            aB2 = aB ^^ twoi
            sB = s + (two ^^ n) * (aB2 - gB ^^ twoi)
            num = 4 * aB2
            den = 1 - sB
        in (num / den) : pis aB gB sB (n + 1)
      puntil f (a:b:xs) = if f a b then b else puntil f (b:xs)
  in puntil (\a b -> abs (a - b) < eps)
     $ pis one (one / sqrt two) zero twoI
 
main :: IO ()
main = do
  -- The last decimal is rounded.
  putStrLn $ toString 1000 $ pi 1000