1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %%%     This is for when it is used in the report      %%%
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %%%         When used as a standalone document         %%%
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 \begin{onlystandalone}
9 \begin{code}
10 {-# LANGUAGE BangPatterns #-}
11 module Mandel where
12 import Complex -- 1.3
13 import PortablePixmap
14 import Control.Parallel
15 import Control.Parallel.Strategies (using)
16 -- import qualified NewStrategies as NS
17 default ()
18 \end{code}
19 \end{onlystandalone}
21         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22 The Mandelbrot set is the set of complex numbers for which the values
23 of the polynomial
24 \begin{equation}
25 \label{poly}
26 f_c(z) = z^2 + c
27 \end{equation}
28 are connected\cite{falconer}. A graphical representation of this set
29 can be rendered by plotting for different points in the complex plane
30 an intensity that is proportional to either the rate of divergence of
31 the above polynomial, or a constant if the polynomial is found to converge.
33 We present a \DPHaskell{} implementation of the mandelbrot set in a bottom up
34 manner that can be decomposed into the following stages :
36 \begin{enumerate}
37 \item   First we generate an infinite list of complex numbers caused by
38         applying the mandelbrot polynomial to a single complex number
39         in the complex plane.
40 \item   Next we walk over the list, determining the rate of divergence of
41         the list of complex numbers.
42 \item   We parallelize the above stages, by mapping the composition
43         of the preceeding functions over the complex plain.
44 \item   Lastly we describe how the complex plain is initialised, and the
45         results are graphically rendered.
46 \end{enumerate}
48         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 \section{Generating the list}
51 \begin{code}
52 mandel::(Num a) => a -> [a]
53 mandel c = infiniteMandel
54            where
55                 infiniteMandel = c : (map (\z -> z*z +c) infiniteMandel)
56 \end{code}
58 The definition of @mandel@ uses the numerically overloaded functions
59 @*@ and @+@; operationally it generates an infinite list of numbers caused by
60 itearting a number (i.e an object that has a type that belongs to class @Num@),
61 with the function @\z -> z*z +c@. For example the evaluation of @mandel 2@''
62 would result in the list of integers :
64 \begin{pseudocode}
65 [2, 6, 38, 1446, 2090918, 4371938082726, 19113842599189892819591078,
66 365338978906606237729724396156395693696687137202086, ^C{Interrupted!}
67 \end{pseudocode}
69 , whereas if the function were applied to the complex number\footnote{complex
70 numbers are defined in Haskells prelude as the algebraic data type
71 @data (RealFloat a) => Complex a = a :+ a@} (1.0 :+ 2.0), then the functions
72 behaviour would be equivalent of the mandelbrot polynomial(\ref{poly}) in
73 which \hbox{$c = 1.0 + i2.0$} :
75 \begin{pseudocode}
76 [(1.0 :+ 2.0), (-2.0 :+ 6.0), (-31.0 :+ -22.0), (478.0 :+ 1366.0),
77 (-1637471.0 :+ 1305898.0), ^C{Interrupted!}
78 \end{pseudocode}
80         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
81 \section{Walking the list}
83 @whenDiverge@'' determines the rate of divergence of a list of complex
84 numbers. If diveregence cannot be proved within a fixed number of
85 iteration of the mandelbrot polynomial, convergence is assumed.
86 This process is encoded in \DPHaskell{} by @take@ing a @limit@ of
87 complex numbers off the infinite stream generated by the mandel function.
88 This {\em finite} list is then traversed, each element of which is checked
89 for divergence.
91 \begin{code}
92 whenDiverge::  Int -> Double -> Complex Double -> Int
94   = walkIt (take limit (mandel c))
95   where
96      walkIt []     = 0                                   -- Converged
97      walkIt (x:xs) | diverge x radius  = 0               -- Diverged
98                    | otherwise           = 1 + walkIt xs -- Keep walking
99 \end{code}
101 \begin{equation}
102 \|x + iy\|_{2} = \sqrt{x^2 + y^2}
103 \end{equation}
105 We assume that divergence occurs if the norm ($\|x\|_{2}$) of a
106 complex number exceeds the radius of the region to be rendered.
107 The predicate @diverge@ encodes the divergence check, where @magnitude@
108 is the prelude function that calculates the euclidean norm
109 ($\|x\|_{2}$) of a complex number.
111 \begin{code}
112 -- VERY IMPORTANT FUNCTION: sits in inner loop
114 diverge::Complex Double -> Double -> Bool
116 \end{code}
118         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119 \section{Parallelising the mandelbrot set}
121 The mandelbrot set can be trivially parallelised by mapping''
122 the @whenDiverge@ function over a complex plain of values.
125 \begin{code}
126 parallelMandel:: [[Complex Double]] -> Int -> Double -> [Int]
128    = concat $129 -- NewStrategies version: 130 -- NS.parListBuffer 50 (NS.seqList id)$
131 --       map (map (whenDiverge limit radius)) mat
133 -- NewStrategies version:
134 --       NS.parListBufferRev 50 (NS.seqList id) $135 -- map (map (whenDiverge limit radius)) mat 137 -- lazyParList version: 138 -- lazyParList 50 139 -- [ let l = map (whenDiverge limit radius) xs 140 -- in seqList l pseq l 141 -- | xs <- mat ] 143 -- lazyParList1 version: 144 parBuffer 70 145 [ let l = map (whenDiverge limit radius) xs 146 in seqList l pseq l 147 | xs <- mat ] 149 -- = lazyParListChunk 100 100$ map (whenDiverge limit radius) mat
150 --   = lazyParMap 512 (whenDiverge limit radius) mat
152 parBuffer :: Int -> [a] -> [a]
153 parBuffer n xs = return xs (start n xs)
154   where
155     return (x:xs) (y:ys) = y par (x : return xs ys)
156     return xs [] = xs
158     start !n [] = []
159     start 0 ys = ys
160     start !n (y:ys) = y par start (n-1) ys
162 -- parListN :: Int -> [a] -> [a]
163 -- parListN 0  xs     = xs
164 -- parListN !n []     = []
165 -- parListN !n (x:xs) = x par parListN (n-1) xs
167 lazyParList :: Int -> [a] -> [a]
168 lazyParList !n xs = go xs (parListN n xs)
169    where
170          go []     _ys    = []
171          go (x:xs) []     = x : xs
172          go (x:xs) (y:ys) = y par (x : go xs ys)
174 lazyParList1 :: Int -> [a] -> [a]
175 lazyParList1 !n xs = go xs (parListN1 n xs [])
176    where
177          go []     _ys    = []
178          go (x:xs) []     = x : xs
179          go (x:xs) (y:ys) = y par (x : go xs ys)
181 -- parMap :: (a -> b) -> [a] -> [b]
182 -- parMap f [] = []
183 -- parMap f (x:xs) = let fx = f x; fxs = parMap f xs in fx par fxs pseq fx:fxs
185 parList :: [a] -> ()
186 parList [] = ()
187 parList (x:xs) = x par parList xs
189 -- parListN version 1: leads to fights as all capabilities try to
190 -- steal the early sparks, and the main thread gets blocked.
191 parListN :: Int -> [a] -> [a]
192 parListN 0  xs     = xs
193 parListN !n []     = []
194 parListN !n (x:xs) = x par parListN (n-1) xs
196 -- like parListN, but starts the sparks in reverse order
197 parListN1 :: Int -> [a] -> [a] -> [a]
198 parListN1 0  xs     ys = parList ys pseq xs
199 parListN1 !n []     ys = parList ys pseq []
200 parListN1 !n (x:xs) ys = parListN1 (n-1) xs (x:ys)
202 seqList :: [a] -> ()
203 seqList [] = ()
204 seqList (x:xs) = x pseq seqList xs
205 --
206 -- parListChunk :: Int -> [a] -> ()
207 -- parListChunk n [] = ()
208 -- parListChunk n xs = let (ys,zs) = splitAt n xs in
209 --                     seqList ys par parListChunk n zs
211 -- parListChunkWHNF :: Int -> [a] -> [a]
212 -- parListChunkWHNF n
213 --   = concat
214 --   . (using parList)
215 --   . map (using seqList)
216 --   . chunk n
218 -- chunk n [] = []
219 -- chunk n xs = as : chunk n bs where (as,bs) = splitAt n xs
220 --
221 -- lazyParList :: Int -> [a] -> [a]
222 -- lazyParList !n xs = go xs (parListN' n xs [])
223 --   where
224 --         go []     _ys    = []
225 --         go (x:xs) []     = x : xs
226 --         go (x:xs) (y:ys) = y par (x : go xs ys)
228 -- lazyParListChunk :: Int -> Int -> [a] -> [a]
229 -- lazyParListChunk !n !size xs = go chunks seqchunks (parListN n seqchunks)
230 --   where
231 --         chunks = chunkList size xs
232 --         seqchunks = map seqList chunks
233 --
234 --         go :: [[a]] -> [()] -> [()] -> [a]
235 --         go []     _ _ys         = []
236 --         go (x:xs) _ []          = concat (x:xs)
237 --         go (x:xs) (y:ys) (z:zs) = z par y pseq (x ++ go xs ys zs)
238 --
239 -- chunkList :: Int -> [a] -> [[a]]
240 -- chunkList !n [] = []
241 -- chunkList !n xs = chunk : chunkList n rest
242 --   where (chunk,rest) = splitAt n xs
243 \end{code}
245 \section{Initialisation of data and graphical rendering.}
247 We render the mandelbrot set by producing a PixMap file that
248 provides a portable representation of a 32bit RGB picture. The format of
249 this file is outside the scope of this paper, but we provide an interface
250 to this representation by defining an abstract data type @Pixmap@ that
251 has the functions @createPixmap@ (signature shown below), and @show@
252 defined (overloaded function from class @Text@).
254 \begin{pseudocode}
255 createPixmap :: Integer ->              -- The width of the Pixmap
256                 Integer ->              -- The height of the Pixmap
257                 Int ->                  -- The depth of of the RGB colour
258                 [(Int,Int,Int)] ->      -- A matrix of the RGB colours
259                 PixMap
260 \end{pseudocode}
262 @mandelset@ controls the evaluation and rendering of the mandelbrot set. The
263 function requires that a viewport'' on the complex plane is defined,
264 such that the region bounded by the viewport will be rendered in a Window''
265 represented by the pixmap.
266 \begin{code}
267 mandelset::Double ->                    -- Minimum X viewport
268            Double ->                    -- Minimum Y viewport
269            Double ->                    -- Maximum X viewport
270            Double ->                    -- maximum Y viewport
271            Integer ->                   -- Window width
272            Integer ->                   -- Window height
273            Int ->                       -- Window depth
274            PixMap                       -- result pixmap
275 mandelset x y x' y' screenX screenY lIMIT
276    = createPixmap screenX screenY lIMIT (map prettyRGB result)
277    where
278 \end{code}
280 @windowToViewport@ is a function that is defined within the scope of the
281 mandelset function (therefore arguments to mandelset will be in scope). The
282 function represents the relationship between the pixels in a window,
283 and points on a complex plain.
285 \begin{code}
286       windowToViewport s t
287            = ((x + (((coerce s) * (x' - x)) / (fromInteger screenX))) :+
288               (y + (((coerce t) * (y' - y)) / (fromInteger screenY))))
290       coerce::Integer -> Double
291       coerce  s   = encodeFloat (toInteger s) 0
292 \end{code}
294 The complex plain is initialised, and the mandelbrot set is calculated.
295 \begin{code}
296       result = parallelMandel
297                   [[windowToViewport s t | s<-[1..screenX]]
298                    | t <- [1..screenY]]
299                   lIMIT
300                   ((max (x'-x) (y'-y)) / 2.0)
302       prettyRGB::Int -> (Int,Int,Int)
303       prettyRGB s = let t = (lIMIT - s) in (s,t,t)
304 \end{code}