Shootout/Pidigits

A Shootout Entry for the pidigits benchmark. It looks like the current entry is affected by this bug in GHC 6.10.1 in which allocating large objects doesn't trigger a GC. Hopefully GHC 6.10.2 will fix this bug and roughly double the speed of the entry.

Proposed Entry

Avoids more multiplications by 0 by ommiting the third element of F all together. A bit faster, and this makes comp1 and comp2 identical again making the code a bit shorter too.

{-# OPTIONS -O2 -optc-O3 #-}
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
-- by Don Stewart, Einar Karttunen, Branimir Maksimovic, 
-- Bertram Felgenhauer, and Darren Smith
--

import System

data F = F !Integer !Integer !Integer

main = loop 10 0 . flip take (str (F 1 0 1) ns) . read . head =<< getArgs

ns = [ F k (4*k+2) (2*k+1) | k <- [1..] ]

loop n s []     = putStrLn $ replicate n ' ' ++ "\t:" ++ show s
loop 0 s xs     = putStrLn ("\t:"++show s) >> loop 10 s xs
loop n s (x:xs) = putStr (show x)          >> loop (n-1) (s+1) xs

flr  x           (F q r t) = (q*x + r) `div` t
comp (F q r t) (F u v x) = F (q*u) (q*v+r*x) (t*x)

str z (x:xs) | y == flr 4 z = y : str (comp (F 10 (-10*y) 1) z) (x:xs)
             | otherwise    =     str (comp z x) xs     where y = flr 3 z

Current Entry

Avoid multiplications by 0 (as seen in C entry). A bit faster.

{-# OPTIONS -O2 -optc-O3 #-}
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
-- by Don Stewart, Einar Karttunen, Branimir Maksimovic and Bertram Felgenhauer
--

import System

data F = F !Integer !Integer !Integer !Integer

main = loop 10 0 . flip take (str (F 1 0 0 1) ns) . read . head =<< getArgs

ns = [ F k (4*k+2) 0 (2*k+1) | k <- [1..] ]

loop n s []     = putStrLn $ replicate n ' ' ++ "\t:" ++ show s
loop 0 s xs     = putStrLn ("\t:"++show s) >> loop 10 s xs
loop n s (x:xs) = putStr (show x)          >> loop (n-1) (s+1) xs

flr  x           (F q r s t) = (q*x + r) `div` (s*x + t)
comp1 (F q r s t) (F u v w x) = F (q*u+r*w) (q*v+r*x) (t*w) (t*x)
comp2 (F q r s t) (F u v w x) = F (q*u) (q*v+r*x) (s*u) (s*v+t*x)

str z (x:xs) | y == flr 4 z = y : str (comp1 (F 10 (-10*y) 0 1) z) (x:xs)  
             | otherwise    =     str (comp2 z x) xs     where y = flr 3 z

Old Entry

Shorter still.

{-# OPTIONS -O2 -optc-O3 #-}
--
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
-- by Don Stewart, Einar Karttunen and Branimir Maksimovic
--

import System

data F = F !Integer !Integer !Integer !Integer

main = loop 10 0 . flip take (str (F 1 0 0 1) ns) . read . head =<< getArgs

ns = [ F k (4*k+2) 0 (2*k+1) | k <- [1..] ]

loop n s []     = putStrLn $ replicate n ' ' ++ "\t:" ++ show s
loop 0 s xs     = putStrLn ("\t:"++show s) >> loop 10 s xs
loop n s (x:xs) = putStr (show x)          >> loop (n-1) (s+1) xs

flr  x           (F q r s t) = (q*x + r) `div` (s*x + t)
comp (F q r s t) (F u v w x) = F (q*u+r*w) (q*v+r*x) (s*u+t*w) (s*v+t*x)

str z (x:xs) | y == flr 4 z = y : str (comp (F 10 (-10*y) 0 1) z) (x:xs)  
             | otherwise    =     str (comp z x) xs     where y = flr 3 z

Older entry

A shorter combination of the Einar's and Branimir's two entries. Same performance. Compile with -O2. At 15 loc, it is by a very long way, the shortest entry, with really good performance. Lesson: arbitrary precsision arithmetic sucks in many other languages.

import System

data LFT = LFT !Integer !Integer !Integer !Integer

floorEx x (LFT q r s t) = (f q * f x + f r) `div` (f s * f x + f t)
    where f = fromInteger

comp (LFT q r s t) (LFT u v w x) = LFT (q*u+r*w) (q*v+r*x) (s*u+t*w) (s*v+t*x)

pi = stream (LFT 1 0 0 1) lfts (floorEx 3) ((==).floorEx 4) prod comp
    where prod z n = comp (LFT 10 (-10*n) 0 1) z
          lfts     = [ LFT k (4*k+2) 0 (2*k+1) | k <- [1..] ]

stream z (x:xs) f g h i
    | g z y     = y : stream (h z y) (x:xs) f g h i 
    | otherwise =     stream (i z x) xs     f g h i where y = f z

main = getArgs >>= loop 10 0 . flip take Main.pi . read . head
  where loop n sum []     = putStrLn $ replicate n ' ' ++ "\t:" ++ show sum
        loop 0 sum xs     = putStrLn ("\t:"++show sum) >> loop 10 sum xs
        loop n sum (x:xs) = putStr (show x)            >> loop (n-1) (sum+1) xs

Original entry

musasabi's entry:

{-# OPTIONS -O2 -optc-O3 #-}
-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
-- contributed by Einar Karttunen
-- adapted from the OCaml version.

import System

floor_ev (q, r, s, t) x = (q*x + r) `div` (s*x + t)
comp (q,r,s,t) (q',r',s',t') = (q*q' + r*s', q*r' + r*t', s*q' + t*s', s*r' + t*t')
next z = floor_ev z 3
safe z n = n == floor_ev z 4
prod z n = comp (10,-10 * n, 0, 1) z
cons z k = let den = 2*k+1 in comp z (fromIntegral k, fromIntegral (2*den), 0, fromIntegral den)

digit :: Int -> (Integer,Integer,Integer,Integer) -> Int -> Int -> Int -> IO ()
digit k z 0 row col = putStrLn (take (10-col) "               "++"\t:"++show (row+col))
digit k z n row col =
  if safe z y
     then if col == 10
         then do let row' = row + 10
                 putStr ("\t:"++show row'++"\n"++show y)
                 digit k (prod z y) (n-1) row' 1
         else putStr (show y) >> digit k (prod z y) (n-1) row (col+1)
     else digit (k+1) (cons z k) n row col
  where y = next z

main = do [n] <- getArgs
          digit 1 (1,0,0,1) (read n) 0 0

Branimir's Idea

An alternative implementation submitted by Branimir Maksimovic

{-# OPTIONS -O2 -optc-O3 #-}
{- original Haskell implementation from spigot.pdf document
 - I've just added printPi and main, also replaced floor and / with div
 - because for some reason div is much faster
 -}
module Main where
import System
main = do [n] <- getArgs
          printPi $ take (read n) Main.pi

printPi digits = printPi' digits 10 0
  where printPi' [] ndigs sum' = do mapM_ (\_ -> putChar ' ') [1..ndigs]
                                    putStr $ "\t:" ++ show sum' ++ "\n"
        printPi' xxs 0 sum' = do putStr $ "\t:" ++ show sum' ++ "\n"
                                 printPi' xxs 10 sum'
        printPi' (x:xs) ndigs sum' = do putStr $ show x
                                        printPi' xs (ndigs-1) (sum'+1)

stream :: (b->c) -> (b->c->Bool) -> (b->c->b) -> (b->a->b) -> b -> [a] -> [c] 
stream next safe prod cons z (x:xs) 
  = if safe z y 
       then y : stream next safe prod cons (prod z y) (x:xs) 
       else stream next safe prod cons (cons z x) xs 
  where y = next z

type LFT = (Integer, Integer, Integer, Integer) 
floorExtr :: LFT -> Integer -> Integer
floorExtr (q,r,s,t) x = ((fromInteger q) * fromInteger x + (fromInteger r)) `div`
                        ((fromInteger s) * fromInteger x + (fromInteger t)) 
unit :: LFT 
unit = (1,0,0,1) 
comp :: LFT -> LFT -> LFT 
comp (q,r,s,t) (u,v,w,x) = (q*u+r*w,q*v+r*x,s*u+t*w,s*v+t*x)

pi = stream next safe prod cons init lfts 
  where
        init = unit 
        lfts = [(k, 4*k+2, 0, 2*k+1) | k<-[1..]] 
        next z = floorExtr z 3
        safe z n = (n == floorExtr z 4) 
        prod z n = comp (10, -10*n, 0, 1) z 
        cons z z'  = comp z z'