Mandelbrot Set Fractal in Fortress

One of my favorite embarrassingly parallel problems is the rendering of the Mandelbrot set fractal. It is a simple program that produces very interesting images.

There are many ways to optimize this program, however I'm going to avoid these optimizations and try to create a naive implementation of the problem, just to see to the program looks like in Fortress. Also I'm going to try to take advantage of the parallel for loops that Fortress provides.


The first thing to define is something to represent complex numbers (remember this is a non-optimize version of the problem) . It seems that in the current reference implementation of there's no native support complex number (that I could find). So this gave me the opportunity to implement it.

value object Complex( real:RR64, img:RR64)
   opr+(self,other:Complex) = Complex(real+other.real,img+other.img)
   opr juxtaposition(self,other:Complex) = 
          Complex((real other.real) - (img other.img),
                  (real other.img) + (other.real img))
   toString() = ("Complex(" real ", " img ")")
end

opr |x:Complex| = SQRT(x.real^2 + x.img^2)


complexExp (o : Complex,n : Integral) = do
   if (n = 1) 
      then o
      else (o complexExp(o,n - 1))
   end
end

opr^(o:Complex,n:Integral) = 
   if (n = 2) 
     then (o o)
     else complexExp(o,n)
   end 

Note that this implementation is not complete, only the required elements for the Mandelbrot set are implemented. It is interesting that multiplication in Fortress doesn't use the '*' operation but it uses juxtaposition of elements to be multiplied. That is, instead of saying (x*y) you say (x y).

Then we need something that helps us to convert from screen coordinates to real coordinates. To do this I created the following function:

lFunc(x1 : RR64,y1 : RR64, x2 : RR64, y2 :RR64) = do
   m = (y2 - y1) / (x2 - x1)
   b = y1 - (m x1)
   fn(x) => (m x) + b
end

Note that this function returns another function that preforms the conversion.

Now the following code shows the implementation of the Mandelbrot set algorithm for one row in the image.

lineMandelbrot[\W\](startIndex : ZZ32, endIndex : ZZ32, 
                    lineData :Array[\ZZ32,ZZ32\], 
                    f :RR64 -> RR64, y : RR64) = do

   for i <- startIndex:endIndex do
      p0 = Complex( f(i) , y)
      p : Complex := p0 
      j : ZZ32 := 0
      maxIteration = 255
      while ( |p| < 2.0  AND j < maxIteration ) do
         p := p^2 + p0
         j := j+1
      end      
      lineData[i] := j
   end
end

Note that I use a parallel for loop for the external for statement in order to say that each point can be calculated independently in a separate thread.

This code looks very nice when formated with Fortify:



Finally the following code show the main program. For the graphic generation, a PPM file was used since is the easiest image format to generate!.

run(a:String...) = 
  do
     imageWidth = 100
     imageHeight = 100
     startY = -1.0
     endY = 1.0
     startX = -1.0
     endX = 1.0

     image =  array[\ZZ32\](imageWidth)

     fout: BufferedWriter = outFileOpen "output.ppm"
     outFileWrite(fout,"P3\n")
     outFileWrite(fout,"" imageWidth " " imageHeight "\n")
     outFileWrite(fout,"255\n")

     
     startIY = 1
     endIY = imageHeight

     fY = lFunc(startIY,startY,endIY,endY)
     fX = lFunc(0,startX,imageWidth - 1,endX)

     for iY <- seq(startIY#endIY) do         
         print "Line: " iY "\n"
         lineMandelbrot(0,imageWidth - 1,image,fX,fY(iY))    

         for i <- seq(0#imageWidth) do
            outFileWrite (fout,"0 0 " image[i]  "   ")
         end
         outFileWrite(fout,"\n")
     end     
     outFileClose(fout)
     ()
end

Since the reference implementation focus is not performance, I think is not fair to publish benchmarks or something like that. However I noticed that by playing around with the NumFortressThreads environment variable (found by looking at the source code) I got different execution times in my dual core machine .