Ticket #1890: B.hs

File B.hs, 2.3 KB (added by dons, 7 years ago)

mandelbrot.hs

Line 
1{-# OPTIONS -fexcess-precision #-}
2--
3-- The Computer Language Shootout
4-- http://shootout.alioth.debian.org/
5--
6-- Contributed by Spencer Janssen, Trevor McCort, Christophe Poucet and Don Stewart
7--
8-- Must be compiled with the -fexcess-precision flag as a pragma. GHC
9-- currently doesn't recognise the -fexcess-precision flag on the command
10-- line (!).
11--
12-- The following flags are suggested when compiling:
13--
14--      -O -fglasgow-exts -optc-march=pentium4
15--      -fbang-patterns -funbox-strict-fields -optc-O2 -optc-mfpmath=sse -optc-msse2
16--
17
18import System
19import System.IO
20import Foreign
21import Foreign.Marshal.Array
22
23main = do
24    w <- getArgs >>= readIO . head
25    let n      = w `div` 8
26        m  = 2 / fromIntegral w
27    putStrLn ("P4\n"++show w++" "++show w)
28    p <- mallocArray0 n
29    unfold n (next_x w m n) p (T 1 0 0 (-1))
30
31------------------------------------------------------------------------
32--
33-- compiled quite differently with ghc 6.8
34--
35-- This function is very sensitive to inlining
36--
37
38unfold :: Int -> (T -> Maybe (Word8,T)) -> Ptr Word8 -> T -> IO ()
39unfold !i !f !ptr !x0 = loop x0
40  where
41    loop !x = go ptr 0 x
42
43    go !p !n !x = case f x of
44        Just (w,y) | n /= i -> poke p w >> go (p `plusPtr` 1) (n+1) y
45        Nothing             -> hPutBuf stdout ptr i
46        _                   -> hPutBuf stdout ptr i >> loop x
47{-# NOINLINE unfold #-}
48
49------------------------------------------------------------------------
50
51--
52-- These are all compiled the same:
53--
54
55data T = T !Int !Int !Int !Double
56
57next_x :: Int -> Double -> Int -> T -> Maybe (Word8, T)
58next_x !w !iw !bw (T bx x y ci)
59    | y  == w   = Nothing
60    | bx == bw  = Just (loop_x w x 8 iw ci 0, T 1 0    (y+1)   (iw+ci))
61    | otherwise = Just (loop_x w x 8 iw ci 0, T (bx+1) (x+8) y ci)
62
63loop_x :: Int -> Int -> Int -> Double -> Double -> Word8 -> Word8
64loop_x !w !x !n !iw !ci !b
65    | x < w = if n == 0
66                    then b
67                    else loop_x w (x+1) (n-1) iw ci (b+b+v)
68    | otherwise = b `shiftL` n
69  where
70    v = fractal 0 0 (fromIntegral x * iw - 1.5) ci 50
71
72fractal :: Double -> Double -> Double -> Double -> Int -> Word8
73fractal !r !i !cr !ci !k
74    | r2 + i2 > 4 = 0
75    | k == 0      = 1
76    | otherwise   = fractal (r2-i2+cr) ((r+r)*i+ci) cr ci (k-1)
77  where
78    (!r2,!i2) = (r*r,i*i)