Personal tools

Shootout/Mandelbrot

From HaskellWiki

< Shootout(Difference between revisions)
Jump to: navigation, search
(Proposed entry)
(Proposed entry)
Line 3: Line 3:
   
 
== Proposed entry ==
 
== Proposed entry ==
  +
  +
Note the explicit disabling of -optc-O3, and turning on sse stuff. Makes
  +
a huge difference on a Pentium M, and a P4.
  +
  +
[https://alioth.debian.org/tracker/index.php?func=detail&aid=304450&group_id=30402&atid=411646 Submitted].
  +
  +
<haskell>
  +
{-# OPTIONS -O -fglasgow-exts -fbang-patterns -funbox-strict-fields -fexcess-precision -optc-O2 -optc-march=pentium4 -optc-mfpmath=sse -optc-msse2 #-}
  +
--
  +
-- The Computer Language Shootout
  +
-- http://shootout.alioth.debian.org/
  +
--
  +
-- Contributed by Spencer Janssen, Trevor McCort, Christophe Poucet and Don Stewart
  +
--
  +
  +
import System
  +
import System.IO
  +
import Foreign
  +
import Foreign.Marshal.Array
  +
  +
main = do
  +
w <- getArgs >>= readIO . head
  +
let n = w `div` 8
  +
m = 2 / fromIntegral w
  +
putStrLn ("P4\n"++show w++" "++show w)
  +
p <- mallocArray0 n
  +
unfold n (next_x w m n) p (T 1 0 0 (-1))
  +
  +
unfold :: Int -> (T -> Maybe (Word8,T)) -> Ptr Word8 -> T -> IO ()
  +
unfold !i !f !ptr !x0 = loop x0
  +
where
  +
loop !x = go ptr 0 x
  +
  +
go !p !n !x = case f x of
  +
Just (w,y) | n /= i -> poke p w >> go (p `plusPtr` 1) (n+1) y
  +
Nothing -> hPutBuf stdout ptr i
  +
_ -> hPutBuf stdout ptr i >> loop x
  +
{-# NOINLINE unfold #-}
  +
  +
data T = T !Int !Int !Int !Double
  +
  +
next_x !w !iw !bw (T bx x y ci)
  +
| y == w = Nothing
  +
| bx == bw = Just (loop_x w x 8 iw ci 0, T 1 0 (y+1) (iw+ci))
  +
| otherwise = Just (loop_x w x 8 iw ci 0, T (bx+1) (x+8) y ci)
  +
  +
loop_x !w !x !n !iw !ci !b
  +
| x < w = if n == 0
  +
then b
  +
else loop_x w (x+1) (n-1) iw ci (b+b+v)
  +
| otherwise = b `shiftL` n
  +
where
  +
v = fractal 0 0 (fromIntegral x * iw - 1.5) ci 50
  +
  +
fractal :: Double -> Double -> Double -> Double -> Int -> Word8
  +
fractal !r !i !cr !ci !k
  +
| r2 + i2 > 4 = 0
  +
| k == 0 = 1
  +
| otherwise = fractal (r2-i2+cr) ((r+r)*i+ci) cr ci (k-1)
  +
where
  +
(!r2,!i2) = (r*r,i*i)
  +
</haskell>
  +
  +
== Old entry ==
   
 
[https://alioth.debian.org/tracker/index.php?func=detail&aid=304447&group_id=30402&atid=411646 Submitted]
 
[https://alioth.debian.org/tracker/index.php?func=detail&aid=304447&group_id=30402&atid=411646 Submitted]
  +
  +
Faster. But still, but gcc is getting in the way!
   
 
<haskell>
 
<haskell>

Revision as of 02:49, 8 February 2007

A Shootout Entry for the mandelbrot benchmark.

Contents

1 Proposed entry

Note the explicit disabling of -optc-O3, and turning on sse stuff. Makes a huge difference on a Pentium M, and a P4.

Submitted.

{-# OPTIONS -O -fglasgow-exts -fbang-patterns -funbox-strict-fields -fexcess-precision -optc-O2 -optc-march=pentium4 -optc-mfpmath=sse -optc-msse2 #-}
--
-- The Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Contributed by Spencer Janssen, Trevor McCort, Christophe Poucet and Don Stewart
--
 
import System
import System.IO
import Foreign
import Foreign.Marshal.Array
 
main = do
    w <- getArgs >>= readIO . head
    let n      = w `div` 8
        m  = 2 / fromIntegral w
    putStrLn ("P4\n"++show w++" "++show w)
    p <- mallocArray0 n
    unfold n (next_x w m n) p (T 1 0 0 (-1))
 
unfold :: Int -> (T -> Maybe (Word8,T)) -> Ptr Word8 -> T -> IO ()
unfold !i !f !ptr !x0 = loop x0
  where
    loop !x = go ptr 0 x
 
    go !p !n !x = case f x of
        Just (w,y) | n /= i -> poke p w >> go (p `plusPtr` 1) (n+1) y
        Nothing             -> hPutBuf stdout ptr i
        _                   -> hPutBuf stdout ptr i >> loop x
{-# NOINLINE unfold #-}
 
data T = T !Int !Int !Int !Double
 
next_x !w !iw !bw (T bx x y ci)
    | y  == w   = Nothing
    | bx == bw  = Just (loop_x w x 8 iw ci 0, T 1 0    (y+1)   (iw+ci))
    | otherwise = Just (loop_x w x 8 iw ci 0, T (bx+1) (x+8) y ci)
 
loop_x !w !x !n !iw !ci !b
    | x < w = if n == 0
                    then b
                    else loop_x w (x+1) (n-1) iw ci (b+b+v)
    | otherwise = b `shiftL` n
  where
    v = fractal 0 0 (fromIntegral x * iw - 1.5) ci 50
 
fractal :: Double -> Double -> Double -> Double -> Int -> Word8
fractal !r !i !cr !ci !k
    | r2 + i2 > 4 = 0
    | k == 0      = 1
    | otherwise   = fractal (r2-i2+cr) ((r+r)*i+ci) cr ci (k-1)
  where
    (!r2,!i2) = (r*r,i*i)

2 Old entry

Submitted

Faster. But still, but gcc is getting in the way!

{-# OPTIONS -fbang-patterns -funbox-strict-fields #-}
--
-- The Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Contributed by Trevor McCort, Spencer Janssen and Don Stewart
-- For best results compile with:
--
-- ghc -O3 -fglasgow-exts -optc-ffast-math -optc-O3 -optc-march=pentium4 -fexcess-precision 
--
 
import System
import Foreign
import qualified Data.ByteString as B
 
main = do
    w <- getArgs >>= readIO . head
    let n      = w `div` 8
        loop_y = B.unfoldrN n (next_x w (2/fromIntegral w) n)
 
        unfold x = case loop_y x of
                    (s, Nothing) -> B.putStr s
                    (s, Just x)  -> B.putStr s >> unfold x
 
    putStrLn ("P4\n"++show w++" "++show w)
    unfold (T 1 0 0 (-1))
 
data T = T !Int !Int !Int !Double
 
next_x !w !iw !bw (T bx x y ci)
    | y  == w   = Nothing
    | bx == bw  = Just (loop_x w x 8 iw ci 0, T 1 0    (y+1)   (iw+ci))
    | otherwise = Just (loop_x w x 8 iw ci 0, T (bx+1) (x+8) y ci)
 
loop_x !w !x !n !iw !ci !b
    | x < w = if n == 0
                    then b
                    else loop_x w (x+1) (n-1) iw ci (b+b+v)
    | otherwise = b `shiftL` n
  where
    v = fractal 0 0 (fromIntegral x * iw - 1.5) ci 50
 
fractal :: Double -> Double -> Double -> Double -> Int -> Word8
fractal !r !i !cr !ci !k
    | r2 + i2 > 4 = 0
    | k == 0      = 1
    | otherwise   = fractal (r2-i2+cr) ((r+r)*i+ci) cr ci (k-1)
  where
    (!r2,!i2) = (r*r,i*i)

3 Current entry

Submitted. Strangely, seems to run slowly on the shootout box.

{-# OPTIONS -fbang-patterns #-}
--
-- The Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Contributed by Trevor McCort, Spencer Janssen and Don Stewart
-- For best results compile with:
--
-- ghc -O3 -fglasgow-exts -optc-ffast-math -optc-O3 -optc-march=pentium4 -fexcess-precision 
--
 
import System
import Foreign
import qualified Data.ByteString.Lazy as B
 
main = do
    (!w) <- getArgs >>= readIO . head
 
    let sh = show $ fromEnum w
        !bw = ceiling (w / 8) :: Int
        !iw = 2/w
 
        gb !ci !x !b !n
            | x == w    = b `shiftL` n
            | n == 0    = b
            | otherwise = gb ci (x+1) (b+b+(lp 0 0 (x * iw - 1.5) ci 50)) (n-1)
 
        ms (bx, x, y, ci)
            | y  == w   = Nothing
            | bx == bw  = Just (gb ci x 0 8,(1,0,y+1, iw+ci))
            | otherwise = Just (gb ci x 0 8,(bx+1,x+8,y,ci))
 
    putStrLn ("P4\n"++sh++" "++sh)
    B.putStr (B.unfoldr ms (1, 0, 0, (-1)))
 
lp :: Double -> Double -> Double -> Double -> Int -> Word8
lp !r !i !cr !ci !k
    | r2 + i2 > 4  = 0
    | k ==  0      = 1
    | otherwise    = lp (r2-i2+cr) ((r+r)*i+ci) cr ci (k-1)
  where
    (!r2,!i2) = (r*r,i*i)
{-# INLINE lp #-}

4 Old entry

-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
-- Based on version by Don Stewart
-- Contributed by Trevor McCort
 
import System
import Data.Bits
import Data.Word
import GHC.Base
 
main = do
    w <- getArgs >>= readIO . head
 
    let ch = chr.fromIntegral
        sh = show $ fromEnum w
        (bw::Int) = ceiling $ w / 8
 
        gb x ci b n
            | x == w    = ch $ b `shiftL` n
            | n == 0    = ch b
            | otherwise = gb (x+1) ci (b+b+(lp 0.0 0.0 50 cr ci)) (n-1)
            where cr = x * 2.0 / w - 1.5
 
        ms bx x y ci
            | y == w    = []
            | bx == bw  = gb x ci 0 8 : ms 1 0 (y+1) ((y+1) * 2.0 / w - 1.0)
            | otherwise = gb x ci 0 8 : ms (bx+1) (x+8) y ci
 
    putStrLn ("P4\n"++sh++" "++sh)
    mapM_ putChar $ ms 1 0 0 (-1.0)
 
lp r i k cr ci | r2 + i2 > (4.0 :: Double) = 0 :: Word32
               | k == (0 :: Word32)        = 1
               | otherwise                 = lp (r2-i2+cr) ((r+r)*i+ci) (k-1) cr ci
    where r2 = r*r ; i2 = i*i

5 Current Entry

Shortest entry in any language.

As with all programs using doubles, compile with -fexcess-precision for big speedups.

-- The Great Computer Language Shootout
-- http://shootout.alioth.debian.org/
-- Based on the SML version, written by Matthias Blume.
-- Implemented in Haskell by Don Stewart
--
import System; import Data.Bits; import Data.Word; import GHC.Base
 
main = do (w::Word32) <- getArgs >>= readIO . head
          putStrLn ("P4\n"++show w++" "++show w) >> yl 0 w w
 
yl y h w = if y < h then xl 0 y 0 8 h w else return ()
 
xl x y b n h w
    | x == w    = putChar (unsafeChr $ b `shiftL` n) >> yl (y+1) h w
    | otherwise = do
        (b',n') <- if n == 0 then putChar (chr b) >> return (0,8) else return (b,n)
        xl (x+1) y (b'+b'+ fromEnum (p x y w h)) (n'-1) h w
 
p (x::Word32) y w h = lp 0.0 0.0 50 (f x * 2.0 / f w - 1.5) (f y * 2.0 / f h - 1.0) 
    where f = fromIntegral
 
lp r i k cr ci | r2 + i2 > (4.0 :: Double) = 0 :: Word32
               | k == (0 :: Word32)        = 1
               | otherwise                 = lp (r2-i2+cr) ((r+r)*i+ci) (k-1) cr ci
    where r2 = r*r ; i2 = i*i

6 Current Entry

The old entry below is 1.2x slower than this version.

This is a translation of the fast SML version. Additionally, we get some good gains by using Word32. (I wonder if this will apply elsewhere?) The -optc-O2 helps as well (another thing to keep in mind for other entries).

{-# OPTIONS -O2 -optc-O2 #-}
--
-- Based on the SML version, written by Matthias Blume.
-- Implemented in Haskell by Don Stewart
--
 
import System
import Data.Bits
import Data.Word
import GHC.Base
 
main = do w <- getArgs >>= return . read . head
          putStrLn $ "P4\n" ++ show w ++ " " ++ show w
          yl 0 w w
 
yl y h w = if y < h then xl 0 y 0 8 h w else return ()
 
xl x y b n h w
    | x == w    = putChar (unsafeChr $ b `shiftL` n) >> yl (y+1) h w
    | otherwise = do
        (b',n') <- if n == 0 then putChar (chr b) >> return (0,8) else return (b,n)
        xl (x+1) y (b'+b'+ fromEnum (p x y w h)) (n'-1) h w
 
p :: Word32 -> Word32 -> Word32 -> Word32 -> Word32
p x y w h = lp 0.0 0.0 50 (f x * 2.0 / f w - 1.5) (f y * 2.0 / f h - 1.0)
    where f = fromIntegral
 
lp r i k cr ci | r2 + i2 > (4.0 :: Double) = 0 :: Word32
               | k == (0 :: Word32)        = 1
               | otherwise                 = lp (r2-i2+cr) ((r+r)*i+ci) (k-1) cr ci
    where (r2,i2) = (r*r, i*i)

7 Original entry

Quite good, though all the lists seem a bit worrying. Also, is putStr legal in this entry?

-- contributed by Greg Buchholz
-- modified by Alson Kemp
-- improvements by Jean-Philippe Bernardy
 
-- compile:  ghc -O2 -o mandelbrot mandelbrot.hs
-- run: mandelbrot 600 >mandel.pbm
 
import Complex
import System(getArgs)
import Char(chr)
import System.IO
 
limit  = 4.0::Double
 
iter   = 50::Int
 
main = do [arg] <- getArgs
          let width = read arg
          --AK:optional;prevent newline mangle on PC
          hSetBinaryMode stdout True
          putStr $ "P4\n" ++ show width ++ " " ++ show width ++ "\n"
          mapM_ putStr $ map (makePBM 0 0) $ fractal (points width width)
 
points :: Int -> Int -> [[Complex Double]]
points width height = [[(2.0*x/w - 1.5) :+ (2.0*y/h - 1) | x<-[0..w-1]] | y<-[0..h-1]]
    where w = fromIntegral width
          h = fromIntegral height
 
fractal :: [[Complex Double]] -> [[Int]]
fractal = map $ map $ fractal' (0.0 :+ 0.0) iter
 
-- magnitude is sloooooowwwwww, so hand code abs^2
fractal' :: Complex Double -> Int -> Complex Double -> Int
fractal' z i c | (realPart z')*(realPart z') + (imagPart z')*(imagPart z') > limit = 0
               | (i == 1) = 1
               | otherwise = fractal' z' (i-1) c
    where z' = z*z+c
 
 
makePBM :: Int -> Int -> [Int] -> [Char]
makePBM i acc []     = chr (acc * 2^(8-i)) : []
makePBM i acc (x:xs) | i==8      = chr acc : makePBM 0 0 (x:xs)
                     | otherwise = makePBM (i+1) (acc*2 + x) xs