Fractals

The Obligatory Fractal Example

Doug Stein, Wolfram Research, Inc.
22 February 1993

This Mathematica Notebook compares the following generators of fractals of iterated Sin and Cos functions on the complex plane:

* a pure Mathematica implementation
* an implementation that uses aMathLink template-based C program
* a parallelized MathLink implementation running on several processors

Examples

One calls the external template-based MathLink program "fractals" differently depending upon whether it is to be parallelized or not. This is because template functions are normally invoked by using LinkWrite to send a CallPacket to the external program; Mathematica then immediately performs a LinkRead on that link and "blocks" until the template program returns. In order to successfully parallelize the external computations, one must execute all of the LinkWrites before the first LinkRead is executed.  The code to do this is found in the implementation section at the bottom of this Notebook.

How to use exactly one external MathLink program

In this example and all those following, we'll use the LinkMode->Connect option and assume that a copy of the fractal engine is running and has opened a MathLink connection in Listen mode. By the way, the default LinkMode for Macintoshes is Launch, which is currently usable only if the template program is to be run on the same CPU as Mathematica.

fractalLinks = OpenFractalLinks[1, LinkMode -> Connect]

Note now the functions that are defined:

LinkPatterns[
      fractalLinks[[1]]]

In addition, the template program defines the following functions in Mathematica when it is Installed by OpenFractalLinks:

Information["SinFractalWrite", LongForm -> True]

Information["CosFractalWrite", LongForm -> True]

SinFractalWrite and CosFractalWrite manually write a CallPacket to the appropriate external MathLink function. Now here is one of the computations that will be used later in the timing tests:

fractalPlot[
      sinFractalC[
        fractalLinks[[1]], {{-7, 7, 
            150}, {-7, 7, 150}}], {0, 50}]

CloseFractalLinks[fractalLinks]

How to use one or more external MathLink programs

fractalLinks = OpenFractalLinks[3, LinkMode -> Connect]

Here we have to call ParallelFractal and SinFractalWrite; these conspire to execute all the LinkWrites first and subsequently execute the corresponding LinkReads. This is somewhat less clear than directly calling sinFractalC, but it does allow parallelism.

fractalPlot[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{0.4, 0.42, 250}, {1.69, 1.71, 250}}], {4, 
        20}]

CloseFractalLinks[fractalLinks]

Timing Tests

In this section I do not use fractalPlot; this is done so timing comparisons do not include PostScript generation and rendering which is constant for all the test cases. Note also that whereas Timing gives Mathematica execution time, AbsoluteTime[] gives wall-clock time.

<<LATE-BREAKING NOTE --- Ignore "Timing" since there seems to be a subtle bug in the accumulation of kernel time when an external MathLink process is very time-consuming; the value of AbsoluteTime[] - startTime is meaningful, however.>>

Native Mathematica Code

startTime = AbsoluteTime[];

Timing[
      sinFractalPureM[{{-7, 7, 150}, {-7, 7, 150}}]; ]

AbsoluteTime[] - startTime

Note: I won't even do the next one because it takes a very, very long time.

startTime = AbsoluteTime[];

Timing[
      sinFractalPureM[{{0.4, 0.42, 250}, {1.69, 1.71, 250}}]; ]

AbsoluteTime[] - startTime

One fractal process on a Macintosh IIci

One copy of fractals881 is running on the same IIci as the Mathematica Front End and Kernel.

fractalLinks = OpenFractalLinks[1, LinkMode -> Connect]

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{-7, 7, 150}, {-7, 7, 150}}]; ]

AbsoluteTime[] - startTime

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{0.4, 0.42, 250}, {1.69, 1.71, 250}}]; ]

AbsoluteTime[] - startTime

CloseFractalLinks[fractalLinks];

2 fractal processes on a IIci

In this test, two copies of fractal881 are run on the same processor as Mathematica's Front End and Kernel.

fractalLinks = OpenFractalLinks[2, LinkMode -> Connect]

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{-7, 7, 150}, {-7, 7, 150}}]; ]

AbsoluteTime[] - startTime

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{0.4, 0.42, 250}, {1.69, 1.71, 250}}]; ]

AbsoluteTime[] - startTime

CloseFractalLinks[fractalLinks];

Three fractal process (1 on a Macintosh IIci and 2 on Rocket 33s)

Here we have three copies of fractals881 running; one is running on each of the Radius Rockets. An additional one is running on the host IIci.

fractalLinks = OpenFractalLinks[3, LinkMode -> Connect]

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{-7, 7, 150}, {-7, 7, 150}}]; ]

AbsoluteTime[] - startTime

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{0.4, 0.42, 250}, {1.69, 1.71, 250}}]; ]

AbsoluteTime[] - startTime

CloseFractalLinks[fractalLinks];

Two fractal processes on Rocket 33s

Here, we're only running 2 copies of fractals881; one copy is running on each of the Radius Rockets.

fractalLinks = OpenFractalLinks[2, LinkMode -> Connect]

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{-7, 7, 150}, {-7, 7, 150}}]; ]

AbsoluteTime[] - startTime

startTime = AbsoluteTime[];

Timing[
      ParallelFractal[fractalLinks, 
        SinFractalWrite, {{0.4, 0.42, 250}, {1.69, 1.71, 250}}]; ]

AbsoluteTime[] - startTime

CloseFractalLinks[fractalLinks];

CosFractal - Another Example

Here are some more pretty pictures; these were done with direct calls to cosFractalC. You can use the parallelized approach above by calling ParallelFractal with a parameter of CosFractalWrite instead of SinFractalWrite.

fractalPlot[
      cosFractalC[{{-7, 7, 500}, {-7, 7, 500}}], {0, 
        20}];

fractalPlot[
      cosFractalC[{{-2.0, -1.75, 500}, {1.3, 1.6, 500}}], {2, 
        20}];

Implementation

Pure Mathematica Implementation

This section defines sinFractalPureM and cosFractalPureM. These are implementations of iterated Sin and Cos fractals on the complex plain written purely in Mathematica.

sinFractalPureM

numberOfSinsToEscape[z_] := 
      Block[{escapeCount, nz = N[z]}, 
        For[escapeCount = 0, Abs[Im[nz]] < 10000 && escapeCount < 100, 
          nz = Sin[nz]; ++escapeCount]; escapeCount]

sinFractalPureM[{{ReMin_, ReMax_, ReSteps_}, {ImMin_, ImMax_, 
            ImSteps_}}] := {fractalObj[{{ReMin, ReMax}, {ImMin, ImMax}}, 
          Table[numberOfSinsToEscape[x + y I], {y, ImMin, 
              ImMax, ((ImMax - ImMin)/ImSteps)}, {x, ReMin, 
              ReMax, ((ReMax - ReMin)/ReSteps)}]]}

cosFractalPureM

numberOfCosToEscape[z_] := 
      Block[{escapeCount, nz = N[z]}, 
        For[escapeCount = 0, Abs[Im[nz]] < 10000 && escapeCount < 100, 
          nz = Cos[nz]; ++escapeCount]; escapeCount]

cosFractalPureM[{{ReMin_, ReMax_, ReSteps_}, {ImMin_, ImMax_, 
            ImSteps_}}] := {fractalObj[{{ReMin, ReMax}, {ImMin, ImMax}}, 
          Table[numberOfCosToEscape[x + y I], {y, ImMin, 
              ImMax, ((ImMax - ImMin)/ImSteps)}, {x, ReMin, 
              ReMax, ((ReMax - ReMin)/ReSteps)}]]}

MathLink and C implementation

This section uses a template-based MathLink program that implements sinFractalC and cosFractalC. The included MathLink template file fractals.tm and C file fractals.c will compile without changes on any system with an ANSI C compiler and a copy of mprep (supplied with Mathematica). I've also included a compiled Macintosh version with this Notebook that requires a 68020 or higher CPU along with a 68881 or higher FPU (a 68040 machine satisfies both criteria). Unfortunately, this version of the fractal program doesn't perform the machine-dependent checks for these so it will crash on a machine that doesn't fit the bill. You can, of course, recompile fractals with whatever compiler settings you desire.

OpenFractalLinks

OpenFractalLinks[maxLinks, opts] will Install up to maxLinks MathLink template programs and return the list of links. It will pass opts to Install so the latter can determine how to open the links. It is simply a convenient wrapper around the Install function, which you could call directly.

OpenFractalLinks[maxLinks_Integer, opts___Rule] := 
      Module[{fractalLinks = {}, aLink, linkcount = maxLinks}, 
        While[linkcount-- > 
              0 && (aLink = Install["fractals", opts]) =!= $Failed, 
          AppendTo[fractalLinks, aLink]]; fractalLinks]

ParallelFractal

ParallelFractal uses a list of links to theMathLink template programs to evaluate the function f on the complex domain {{ReMin, ReMax},{ImMin,ImMax}} with an escape parameter of steps. It divides the problem up among the processes as equally as possible. My definition of "equal" is based on area in the complex plain, not execution time.

ParallelFractal[fractalLinks_List, 
        f_, {{ReMin_, ReMax_, ReSteps_}, {ImMin_, ImMax_, ImSteps_}}] := 
      Module[{numLinks = Length[fractalLinks], useTheseLinks, ranges, itvar, 
          pieceSteps, pieces, ans}, 
        If[ImSteps ≤ numLinks, 
          useTheseLinks = Take[fractalLinks, ImSteps], 
          useTheseLinks = fractalLinks]; 
        pieceSteps = Ceiling[(ImSteps/numLinks)]; 
        ranges = 
          Partition[
            Table[itvar, {itvar, ImMin, ImMax, ((ImMax - ImMin)/numLinks)}], 
            2, 1]; MapThread[
          f[#1, {{ReMin, ReMax, ReSteps}, 
                Join[#2, {pieceSteps}]}] &, {useTheseLinks, N[ranges]}]; 
        pieces = LinkRead /@ useTheseLinks; 
        ans = Join @@ (#1[[2]] &) /@ 
              pieces; 
        ans = fractalObj @@ {{{ReMin, ReMax}, {ImMin, ImMax}}, ans}; 
        ans]

CloseFractalLinks

CloseFractalLinks[l_List] := Uninstall /@ l;

fractalPlot

fractalPlot is a wrapper around ListDensityPlot.

fractalPlot[fractal_fractalObj, {minCount_, maxCount_}] := 
      Module[{}, 
        ListDensityPlot[fractal[[2]], 
          PlotRange -> {minCount, maxCount}, 
          MeshRange -> 
            fractal[[1]], 
          Mesh -> False, AspectRatio -> Automatic, 
          ColorFunction -> Hue]]