Speed improvements - optimisation

Symmetry

The Mandelbrot set is symmetric with respect to the x-axis in the plane :

colour(x,y) = colour(x,-y)

its intersection with the x-axis ( real slice of Mandelbrot set ) is an interval :

<-2 ; 1/4>

It can be used to speed up computations ( up to 2-times )^[1]

Bailout test

Instead of checking if magnitude ( radius = abs(z) ) is greater than the escape radius ( ER):

 ${\sqrt {Z_{x}^{2}+Z_{y}^{2}}}<ER$

compute square of ER:

 ER2 = ER*ER

once and check :

 $Z_{x}^{2}+Z_{y}^{2}<ER^{2}$

It gives the same result and is faster.

interior detection

Interior detection method^[2]

time with detection versus without detection is : 0.383s versus 8.692s so it is 23 times faster !!!!

// gives last iterate = escape time
// output 0< i < iMax
 int iterate(double complex C , int iMax)
  {
   int i=0;
   double complex Z= C; // initial value for iteration Z0
   complex double D = 1.0; // derivative with respect to z 
   
   for(i=0;i<iMax;i++)
    { if(cabs(Z)>EscapeRadius) break; // exterior
      if(cabs(D)<eps) break; // interior
      D = 2.0*D*Z;
      Z=Z*Z+C; // complex quadratic polynomial
      
    }
   return i; 
 }

Period detection

Period of Complex_quadratic_map - main article

"When rendering a Mandelbrot or Julia set, the most time-consuming parts of the image are the “black areas”. In these areas, the iterations never escape to infinity, so every pixel must be iterated to the maximum limit. Therefore, much time can be saved by using an algorithm to detect these areas in advance, so that they can be immediately coloured black, rather than rendering them in the normal way, pixel by pixel. FractalNet uses a recursive algorithm to split the image up into blocks, and tests each block to see whether it lies inside a “black area”. In this way, large areas of the image can be quickly coloured black, often saving a lot of rendering time. ... (some) blocks were detected as “black areas” and coloured black immediately, without having to be rendered. (Other) blocks were rendered in the normal way, pixel by pixel." Michael Hogg ^[3]

 // cpp code by Geek3 
// http://commons.wikimedia.org/wiki/File:Mandelbrot_set_rainbow_colors.png
bool outcircle(double center_x, double center_y, double r, double x, double y)
{ // checks if (x,y) is outside the circle around (center_x,center_y) with radius r
        x -= center_x;
        y -= center_y;
        if (x * x + y * y > r * r)
                return(true);
        return(false);

 // skip values we know they are inside
                        if ((outcircle(-0.11, 0.0, 0.63, x0, y0) || x0 > 0.1)
                                && outcircle(-1.0, 0.0, 0.25, x0, y0)
                                && outcircle(-0.125, 0.744, 0.092, x0, y0)
                                && outcircle(-1.308, 0.0, 0.058, x0, y0)
                                && outcircle(0.0, 0.25, 0.35, x0, y0))
                        {
                          // code for iteration
                         }

Cardioid and period-2 checking

One way to improve calculations is to find out beforehand whether the given point lies within the cardioid or in the period-2 bulb. Before passing the complex value through the escape time algorithm, first check if:

(x+1)^{2}+y^{2}<{\frac {1}{16}}

to determine if the point lies within the period-2 bulb and

q=\left(x-{\frac {1}{4}}\right)^{2}+y^{2}

q\left(q+\left(x-{\frac {1}{4}}\right)\right)<{\frac {1}{4}}y^{2}.

to determine if the point lies inside the main cardioid.

another description ^[4]

// http://www.fractalforums.com/new-theories-and-research/quick-(non-iterative)-rejection-filter-for-mandelbrotbuddhabrot-with-benchmark/
         public static void quickRejectionFilter(BlockingCollection<Complex> input, BlockingCollection<Complex> output)
         {
            foreach(var item in input.GetConsumingEnumerable())
            {
                if ((Complex.Abs(1.0 - Complex.Sqrt(Complex.One - (4 * item))) < 1.0)) continue;
                if (((Complex.Abs(item - new Complex(-1, 0))) < 0.25)) continue;
                if ((((item.Real + 1.309) * (item.Real + 1.309)) + item.Imaginary * item.Imaginary) < 0.00345) continue;
                if ((((item.Real + 0.125) * (item.Real + 0.125)) + (item.Imaginary - 0.744) * (item.Imaginary - 0.744)) < 0.0088) continue;
                if ((((item.Real + 0.125) * (item.Real + 0.125)) + (item.Imaginary + 0.744) * (item.Imaginary + 0.744)) < 0.0088) continue;

                //We tried every known quick filter and didn't reject the item, adding it to next queue.
                output.Add(item);
            }
         }

Periodicity checking

Most points inside the Mandelbrot set oscillate within a fixed orbit. There could be anything from ten to tens of thousands of points in between, but if an orbit ever reaches a point where it has been before then it will continually follow this path, will never tend towards infinity and hence is in the Mandelbrot set.

This Mandelbrot processor includes:

periodicity checking
period-2 bulb and cardioid checking

for a great speed up during deep zoom animations with a high maximum iteration value.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Mandelbrot2
{
    public struct MandelbrotData
    {
        private double px;
        private double py;

        private double zx;
        private double zy;

        private long iteration;
        private bool inSet;
        private HowFound found;

        public MandelbrotData(double px, double py)
        {
            this.px = px;
            this.py = py;
            this.zx = 0.0;
            this.zy = 0.0;
            this.iteration = 0L;
            this.inSet = false;
            this.found = HowFound.Not;
        }

        public MandelbrotData(double px, double py,
                              double zx, double zy,
                              long iteration,
                              bool inSet,
                              HowFound found)
        {
            this.px = px;
            this.py = py;
            this.zx = zx;
            this.zy = zy;
            this.iteration = iteration;
            this.inSet = inSet;
            this.found = found;
        }

        public double PX
        {
            get { return this.px; }
        }

        public double PY
        {
            get { return this.py; }
        }

        public double ZX
        {
            get { return this.zx; }
        }

        public double ZY
        {
            get { return this.zy; }
        }

        public long Iteration
        {
            get { return this.iteration; }
        }

        public bool InSet
        {
            get { return this.inSet; }
        }

        public HowFound Found
        {
            get { return this.found; }
        }
    }

    public enum HowFound { Max, Period, Cardioid, Bulb, Not }

    class MandelbrotProcess
    {
        private long maxIteration;
        private double bailout;

        public MandelbrotProcess(long maxIteration, double bailout)
        {
            this.maxIteration = maxIteration;
            this.bailout = bailout;
        }

        public MandelbrotData Process(MandelbrotData data)
        {
            double zx;
            double zy;
            double xx;
            double yy;

            double px = data.PX;
            double py = data.PY;
            yy = py * py;

            #region Cardioid check

            //Cardioid
            double temp = px - 0.25;
            double q = temp * temp + yy;
            double a = q * (q + temp);
            double b = 0.25 * yy;
            if (a < b)
                return new MandelbrotData(px, py, px, py, this.maxIteration, true, HowFound.Cardioid);

            #endregion

            #region Period-2 bulb check

            //Period-2 bulb
            temp = px + 1.0;
            temp = temp * temp + yy;
            if (temp < 0.0625)
                return new MandelbrotData(px, py, px, py, this.maxIteration, true, HowFound.Bulb);

            #endregion

            zx = px;
            zy = py;

            int check = 3;
            int checkCounter = 0;

            int update = 10;
            int updateCounter = 0;

            double hx = 0.0;
            double hy = 0.0;

            for (long i = 1; i <= this.maxIteration; i++)
            {
                //Calculate squares
                xx = zx * zx;
                yy = zy * zy;

                #region Bailout check

                //Check bailout
                if (xx + yy > this.bailout)
                    return new MandelbrotData(px, py, zx, zy, i, false, HowFound.Not);

                #endregion

                //Iterate
                zy = 2.0 * zx * zy + py;
                zx = xx - yy + px;

                #region Periodicity check

                //Check for period
                double xDiff = Math.Abs(zx - hx);
                if (xDiff < this.ZERO)
                {
                    double yDiff = Math.Abs(zy - hy);
                    if (yDiff < this.ZERO)
                        return new MandelbrotData(px, py, zx, zy, i, true, HowFound.Period);
                } //End of on zero tests

                //Update history
                if (check == checkCounter)
                {
                    checkCounter = 0;

                    //Double the value of check
                    if (update == updateCounter)
                    {
                        updateCounter = 0;
                        check *= 2;
                    }
                    updateCounter++;

                    hx = zx;
                    hy = zy;
                } //End of update history
                checkCounter++;

                #endregion
            } //End of iterate for

            #region Max iteration

            return new MandelbrotData(px, py, zx, zy, this.maxIteration, true, HowFound.Max);

            #endregion
        }

        private double ZERO = 1e-17;
    }
}

References

T Myers: Introduction to perturbation ...
interpretation-of-topology-optimization
Derivative Extrapolation From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2023.

[1] How to use symetry of set

[2] A Cheritat wiki: Mandelbrot_set#Interior_detection_methods

[3] FractalNet by Michael Hogg

[4] Fractal forums: quick-(non-iterative)-rejection-filter-for-mandelbrotbuddhabrot-with-benchmark/

[1]

[2]

[3]

[4]