Chudnovskyの公式を用いた円周率の計算

Chudnovsky.hs

import Data.Bits
import Control.Parallel

cutoff = 1000000

rsqrt2 :: Int -> Int -> Integer
rsqrt2 d n = iter 32 $ (floor $ approx * 2 ** 64) `shiftL` (d - 64)
  where
    approx = 1.0 / sqrt (fromIntegral n :: Double)
    mult x y = (x * y) `shiftR` d
    iter p x
      | p >= d = x
      | otherwise = iter (p * 2) ((3 * x - toInteger n * x `mult` x `mult` x) `shiftR` 1)

series :: Integer -> Integer -> (Integer, Integer, Integer)
series l r
  | r - l == 1 =
      let pr = (2 * r - 1) * (6 * r - 5) * (6 * r - 1)
          qr = r * r * r * 10939058860032000
          ar = (if even r then 1 else -1) * (13591409 + 545140134 * r)
      in (pr, qr, ar * pr)
  | cutoff <= r - l =
      let m = (r + l) `div` 2
          s1 = series l m
          s2 = series m r
          (p1, q1, t1) = s1
          (p2, q2, t2) = s2
      in s1 `par` (s2 `pseq` (p1 * p2, q1 * q2, t1 * q2 + p1 * t2))
  | otherwise =
      let m = (r + l) `div` 2
          (p1, q1, t1) = series l m
          (p2, q2, t2) = series m r
      in (p1 * p2, q1 * q2, t1 * q2 + p1 * t2)

main :: IO ()
main = do
  s <- getLine
  let n = read s
      n2 = (+ 1024) $ ceiling $ fromIntegral n * logBase 2 10
      se = series 0 $ (n + 13) `div` 14
      r = rsqrt2 n2 10005
      (_, q, t) = se
  putStrLn $ show $ s `par` (r `pseq` (((4270934400 * r * q) `div` (t + 13591409 * q)) * 10 ^ n) `shiftR` n2)