Personal tools

Shootout/Mandelbrot

From HaskellWiki

< Shootout
Revision as of 04:10, 7 February 2007 by DonStewart (Talk | contribs)

Jump to: navigation, search

A Shootout Entry for the mandelbrot benchmark.

Contents

1 Proposed entry

Submitted.

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

2 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

3 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

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

5 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