[project @ 1996-07-25 21:02:03 by partain]
[nofib.git] / spectral / mandel / Mandel.lhs
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %%%     This is for when it is used in the report      %%%
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4
5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 %%%         When used as a standalone document         %%%
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 \begin{onlystandalone}
9 \begin{code}
10 module Mandel where
11 import Complex -- 1.3
12 import PortablePixmap
13 default ()
14 \end{code}
15 \end{onlystandalone}
16
17         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18 The Mandelbrot set is the set of complex numbers for which the values
19 of the polynomial
20 \begin{equation}
21 \label{poly}
22 f_c(z) = z^2 + c
23 \end{equation}
24 are connected\cite{falconer}. A graphical representation of this set
25 can be rendered by plotting for different points in the complex plane
26 an intensity that is proportional to either the rate of divergence of
27 the above polynomial, or a constant if the polynomial is found to converge.
28
29 We present a \DPHaskell{} implementation of the mandelbrot set in a bottom up
30 manner that can be decomposed into the following stages :
31
32 \begin{enumerate}
33 \item   First we generate an infinite list of complex numbers caused by
34         applying the mandelbrot polynomial to a single complex number
35         in the complex plane.
36 \item   Next we walk over the list, determining the rate of divergence of
37         the list of complex numbers.
38 \item   We parallelize the above stages, by mapping the composition
39         of the preceeding functions over the complex plain.
40 \item   Lastly we describe how the complex plain is initialised, and the
41         results are graphically rendered.
42 \end{enumerate}
43
44         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 \section{Generating the list}
46
47 \begin{code}
48 mandel::(Num a) => a -> [a]
49 mandel c = infiniteMandel
50            where
51                 infiniteMandel = c : (map (\z -> z*z +c) infiniteMandel)
52 \end{code}
53
54 The definition of @mandel@ uses the numerically overloaded functions
55 @*@ and @+@; operationally it generates an infinite list of numbers caused by
56 itearting a number (i.e an object that has a type that belongs to class @Num@),
57 with the function @\z -> z*z +c@. For example the evaluation of ``@mandel 2@''
58 would result in the list of integers :
59
60 \begin{pseudocode}
61 [2, 6, 38, 1446, 2090918, 4371938082726, 19113842599189892819591078,
62 365338978906606237729724396156395693696687137202086, ^C{Interrupted!}
63 \end{pseudocode}
64
65 , whereas if the function were applied to the complex number\footnote{complex
66 numbers are defined in Haskells prelude as the algebraic data type
67 @data (RealFloat a) => Complex a = a :+ a@} (1.0 :+ 2.0), then the functions
68 behaviour would be equivalent of the mandelbrot polynomial(\ref{poly}) in
69 which \hbox{$c = 1.0 + i2.0$} :
70
71 \begin{pseudocode}
72 [(1.0 :+ 2.0), (-2.0 :+ 6.0), (-31.0 :+ -22.0), (478.0 :+ 1366.0),
73 (-1637471.0 :+ 1305898.0), ^C{Interrupted!}
74 \end{pseudocode}
75
76         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77 \section{Walking the list}
78
79 ``@whenDiverge@'' determines the rate of divergence of a list of complex
80 numbers. If diveregence cannot be proved within a fixed number of
81 iteration of the mandelbrot polynomial, convergence is assumed.
82 This process is encoded in \DPHaskell{} by @take@ing a @limit@ of
83 complex numbers off the infinite stream generated by the mandel function.
84 This {\em finite} list is then traversed, each element of which is checked
85 for divergence.
86
87 \begin{code}
88 whenDiverge::  Int -> Double -> Complex Double -> Int
89 whenDiverge limit radius c
90   = walkIt (take limit (mandel c))
91   where
92      walkIt []     = 0                                   -- Converged
93      walkIt (x:xs) | diverge x radius  = 0               -- Diverged
94                    | otherwise           = 1 + walkIt xs -- Keep walking
95 \end{code}
96
97 \begin{equation}
98 \|x + iy\|_{2} = \sqrt{x^2 + y^2}
99 \end{equation}
100
101 We assume that divergence occurs if the norm ($\|x\|_{2}$) of a
102 complex number exceeds the radius of the region to be rendered.
103 The predicate @diverge@ encodes the divergence check, where @magnitude@
104 is the prelude function that calculates the euclidean norm
105 ($\|x\|_{2}$) of a complex number.
106
107 \begin{code}
108 diverge::Complex Double -> Double -> Bool
109 diverge cmplx radius =  magnitude cmplx > radius
110 \end{code}
111
112         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113 \section{Parallelising the mandelbrot set}
114
115 The mandelbrot set can be trivially parallelised by ``mapping''
116 the @whenDiverge@ function over a complex plain of values.
117
118
119 \begin{code}
120 parallelMandel:: [Complex Double] -> Int -> Double -> [Int]
121 parallelMandel mat limit radius
122    = map (whenDiverge limit radius) mat
123 \end{code}
124
125 \section{Initialisation of data and graphical rendering.}
126
127 We render the mandelbrot set by producing a PixMap file that
128 provides a portable representation of a 32bit RGB picture. The format of
129 this file is outside the scope of this paper, but we provide an interface
130 to this representation by defining an abstract data type @Pixmap@ that
131 has the functions @createPixmap@ (signature shown below), and @show@
132 defined (overloaded function from class @Text@).
133
134 \begin{pseudocode}
135 createPixmap :: Integer ->              -- The width of the Pixmap
136                 Integer ->              -- The height of the Pixmap
137                 Int ->                  -- The depth of of the RGB colour
138                 [(Int,Int,Int)] ->      -- A matrix of the RGB colours
139                 PixMap
140 \end{pseudocode}
141
142 @mandelset@ controls the evaluation and rendering of the mandelbrot set. The
143 function requires that a ``viewport'' on the complex plane is defined,
144 such that the region bounded by the viewport will be rendered in a ``Window''
145 represented by the pixmap.
146 \begin{code}
147 mandelset::Double ->                    -- Minimum X viewport
148            Double ->                    -- Minimum Y viewport
149            Double ->                    -- Maximum X viewport
150            Double ->                    -- maximum Y viewport
151            Integer ->                   -- Window width
152            Integer ->                   -- Window height
153            Int ->                       -- Window depth
154            PixMap                       -- result pixmap
155 mandelset x y x' y' screenX screenY lIMIT
156    = createPixmap screenX screenY lIMIT (map prettyRGB result)
157    where
158 \end{code}
159
160 @windowToViewport@ is a function that is defined within the scope of the
161 mandelset function (therefore arguments to mandelset will be in scope). The
162 function represents the relationship between the pixels in a window,
163 and points on a complex plain.
164
165 \begin{code}
166       windowToViewport s t
167            = ((x + (((coerce s) * (x' - x)) / (fromInteger screenX))) :+
168               (y + (((coerce t) * (y' - y)) / (fromInteger screenY))))
169
170       coerce::Integer -> Double
171       coerce  s   = encodeFloat (toInteger s) 0
172 \end{code}
173
174 The complex plain is initialised, and the mandelbrot set is calculated.
175 \begin{code}
176       result = parallelMandel
177                   [windowToViewport s t | t <- [1..screenY] , s<-[1..screenX]]
178                   lIMIT
179                   ((max (x'-x) (y'-y)) / 2.0)
180
181       prettyRGB::Int -> (Int,Int,Int)
182       prettyRGB s = let t = (lIMIT - s) in (s,t,t)
183 \end{code}