Machinの公式を使ってπを1000桁まで計算してみました

円周率 10000 桁を参考にπを1000桁まで求めるプログラムを書いてみました。HaskellだとInteger(桁制限のない整数)とかRational(有理数)があるので意外と楽々でした。級数展開のとこも素直に書けたし。
なおRationalを文字列にする関数は、有理数の値を指定した桁数だけ表示するのものを使わせていただきました。

module Main (main) where

import Data.Ratio

arctan1per :: Integer -> Rational
arctan1per m = atan1 0
  where
    atan1 :: Integer -> Rational
    atan1 n
      | a < permissibleError = 0
      | n `mod` 2 == 0 =  a + atan1 (n + 1)
      | otherwise      = -a + atan1 (n + 1)
      where
        a = let x = 2 * n + 1 in 1 % (x * (m ^ x))

calcPi :: Rational
calcPi = 4 * (4 * arctan1per 5 - arctan1per 239)

significantDigits :: Int
significantDigits = 1000

permissibleError :: Rational
permissibleError = 1 % (10 ^ (significantDigits + 1))

main :: IO ()
main = dump $ calcPi
  where
    showRational :: Int -> Rational -> String
    showRational n r = h ++ ('.' : l)
      where
        dm0 = divMod x y
        x = numerator r
        y = denominator r
        s = take (n + 1) $ iterate (flip divMod y . (10 *) . snd) dm0
        h = show $ fst $ head s
        l = concatMap (show . fst) $ tail s

    dump :: Rational -> IO ()
    dump r = do
      let (xs1, xs2) = splitAt 2 $ showRational significantDigits r
      putStr xs1
      dump' 0 xs2
      putStrLn ""
      where
        dump' _ [] = return ()
        dump' n xs = do
          let (xs1, xs2) = splitAt 4 xs
              sep = if n `mod` 12 == 11 then " \n  " else " "
          putStr $ xs1 ++ sep
          dump' (n + 1) xs2

-- => 3.1415 9265 3589 7932 3846 2643 3832 7950 2884 1971 6939 9375
--      1058 2097 4944 5923 0781 6406 2862 0899 8628 0348 2534 2117
--      0679 8214 8086 5132 8230 6647 0938 4460 9550 5822 3172 5359
--      4081 2848 1117 4502 8410 2701 9385 2110 5559 6446 2294 8954
--      9303 8196 4428 8109 7566 5933 4461 2847 5648 2337 8678 3165
--      2712 0190 9145 6485 6692 3460 3486 1045 4326 6482 1339 3607
--      2602 4914 1273 7245 8700 6606 3155 8817 4881 5209 2096 2829
--      2540 9171 5364 3678 9259 0360 0113 3053 0548 8204 6652 1384
--      1469 5194 1511 6094 3305 7270 3657 5959 1953 0921 8611 7381
--      9326 1179 3105 1185 4807 4462 3799 6274 9567 3518 8575 2724
--      8912 2793 8183 0119 4912 9833 6733 6244 0656 6430 8602 1394
--      9463 9522 4737 1907 0217 9860 9437 0277 0539 2171 7629 3176
--      7523 8467 4818 4676 6940 5132 0005 6812 7145 2635 6082 7785
--      7713 4275 7789 6091 7363 7178 7214 6844 0901 2249 5343 0146
--      5495 8537 1050 7922 7968 9258 9235 4201 9956 1121 2902 1960
--      8640 3441 8159 8136 2977 4771 3099 6051 8707 2113 4999 9998
--      3729 7804 9951 0597 3173 2816 0963 1859 5024 4594 5534 6908
--      3026 4252 2308 2533 4468 5035 2619 3118 8171 0100 0313 7838
--      7528 8658 7533 2083 8142 0617 1776 6914 7303 5982 5349 0428
--      7554 6873 1159 5628 6388 2353 7875 9375 1957 7818 5778 0532
--      1712 2680 6613 0019 2787 6611 1959 0921 6420 1989