Shootout/Mandelbrot

A Shootout Entry for the mandelbrot benchmark.

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 -fexcess-precision -fbang-patterns #-}
--
-- The Computer Language Shootout
-- http://shootout.alioth.debian.org/
--
-- Contributed by Spencer Janssen, Trevor McCort, Christophe Poucet and Don Stewart
--
-- Must be compiled with the -fexcess-precision flag as a pragma. GHC
-- currently doesn't recognise the -fexcess-precision flag on the command
-- line (!).
--
-- The following flags are suggested when compiling:
--
--      -O -fglasgow-exts -optc-march=pentium4
--      -fbang-patterns -funbox-strict-fields -optc-O2 -optc-mfpmath=sse -optc-msse2
--

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))

{-# NOINLINE unfold #-}
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

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 :: Int -> Int -> Int -> Double -> Double -> Word8 -> Word8
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)

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)

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 #-}

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)
    putStr $ 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

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

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)

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