import Control.Monad
import Data.Number.CReal
import GHC.Integer
import Text.Printf

iterations = 52
main = do
  printf "N.%44s%4s%s\n" 
          "Integral part of Nth term" "10^" "=Actual value of Nth term"

  forM_ [0..9] $ \n ->
    printf "%d.%44d%4d%s\n" n
                               (almkvistGiulleraIntegral n)
                               (tenExponent n)
                               (showCReal 50 (almkvistGiullera n))

  printf "\nPi after%d iterations:\n" iterations
  putStrLn $ showCReal 70 $ almkvistGiulleraPi iterations


almkvistGiulleraIntegral n =
  let polynomial  = (532 `timesInteger` n `timesInteger` n) `plusInteger` (126 `timesInteger` n) `plusInteger` 9
      numerator   = 32 `timesInteger` (facInteger (6 `timesInteger` n)) `timesInteger` polynomial
      denominator = 3 `timesInteger` (powInteger (facInteger n) 6)
   in numerator `divInteger` denominator


tenExponent n = 3 `minusInteger` (6 `timesInteger` (1 `plusInteger` n))


almkvistGiullera n = fromInteger (almkvistGiulleraIntegral n) / fromInteger (powInteger 10 (abs (tenExponent n)))


almkvistGiulleraSum n = sum $ map almkvistGiullera [0 .. n]


almkvistGiulleraPi n = sqrt $ 1 / almkvistGiulleraSum n


facInteger n = if n `leInteger` 1 then 1 else n `timesInteger` facInteger (n `minusInteger` 1)


powInteger 1 _ = 1
powInteger _ 0 = 1
powInteger b 1 = b
powInteger b e = b `timesInteger` powInteger b (e `minusInteger` 1)