Help optimising mandelbrot

[PurpleBlu3s]

I am trying to optimise my mandelbrot algorithm using this: http://en.wikipedia.org/wiki/Mandelbrot_set#Optimizations
However I'm stuck on how to implement the cardioid ones (first 2). I have done polar coordinates but it was a while ago and only introductory level. Wikipedia doesn't explain what p and q are, so I'm confused how to implement it. Any help would be greatly appreciated.

Thanks.

[msltoe]

In the link you provided, p and q are simply intermediate variables used for the "if" test. They are defined in the equations as functions of x (real component) and y (imaginary component).

There may be optimized libraries out there for computing the M-set that take advantage of all of the vectorization (parallelization) available on x86 chips. Also, some people on the forum write directly in assembly language to get the best optimization.

[PurpleBlu3s]

In the link you provided, p and q are simply intermediate variables used for the "if" test. They are defined in the equations as functions of x (real component) and y (imaginary component).

There may be optimized libraries out there for computing the M-set that take advantage of all of the vectorization (parallelization) available on x86 chips. Also, some people on the forum write directly in assembly language to get the best optimization.

How then should p and/or q be defined?

[PurpleBlu3s]

Ah, nevermind, I just hadn't read the article carefully.

Didn't realise they were defining p/q in the first equation - I thought that was one of the conditions to be satisfied.

Just quickly tested it - saw a marked improvement. The optimised time was about 67% of the non-optimised time.

[zenzero-2001]

In my opinion, cardioid checking (and testing for the low order buds) isn't a very useful strategy for optimising the calculation of the Mandelbrot fractal. This technique is only good for very low zooms, which only require a low number of iterations and are therefore fast without any optimisation.

More effective techniques for optimising the calculation (that work at higher zoom levels) are periodicity checking (also known as obit detection) and adjacency optimisation. I recommend having a read of Robert Munafo's Encyclopedia of the Mandelbrot Set: http://mrob.com/pub/muency.html . He describes various techniques for adjacency optimisation: Solid guessing (successive refinement), boundary tracing and Mariani / Silver (also known as tesseral). All these techniques will help speed up the calculation inside the set (the black area).

There are two approaches to periodicity checking, one of which is described on Robert's website (the RHOP algorithm, the other being the binary back off algorithm). Both techniques will add overhead to the calculation but will greatly speed up detecting points inside the set. I believe the RHOP algorithm is faster than the BBO algorithm, but adds greater overhead (could be wrong though). I use the latter algorithm in my Mandelbrot program.

I have also found this page with a listing in Java, which you may find useful (includes cardioid / bulb detection and period checking - using the binary back off algorithm):

http://en.wikipedia.org/wiki/User:Simpsons_contributor/periodicity_checking

zenzero-2001

[fractalmind]

I use this approach:

 a = c_re;
    b = c_im;

    for( i=3; i<LIM; i+=4 ) {
        oa = a;
        ob = b;

        x = 2*a*b + c_im;
        a = a*a - b*b + c_re;

        b = 2*a*x + c_im;
        a = a*a - x*x + c_re;

        x = 2*a*b + c_im;
        a = a*a - b*b + c_re;

        b = a*x;
        x *= x;
        a *= a;

        if( a + x > 4.0 ) break;   

        a = a - x + c_re;
        b = b + b + c_em;
    }

    a = oa;
    b = ob;
    i -= 3;

    for( ; i<LIM; i++ ) {
        x = a*b;
        a *= a;
        b *= b;
        if( a + b > 4.0 ) break;   
        a = a - b + c_re;
        b = x + x + c_im;
    }

    plot( c_im, c_re, i );

http://www.azillionmonkeys.com/qed/p5opt.html

By converting this c-code to assembler some additional speed is gained. Also it is better to use bitblitting (like windows BITBLT) instructions to render many pixels at the same time instead of plotting them individually)

[ker2x]

I am trying to optimise my mandelbrot algorithm using this: http://en.wikipedia.org/wiki/Mandelbrot_set#Optimizations
However I'm stuck on how to implement the cardioid ones (first 2). I have done polar coordinates but it was a while ago and only introductory level. Wikipedia doesn't explain what p and q are, so I'm confused how to implement it. Any help would be greatly appreciated.

Thanks.

i had the same problem, please take a look at this thread : http://www.fractalforums.com/programming/need-help-to-convert-(abs(-1-0-sqrt(1-(4*c))-)-to-c/

And my OpenCL implementation of a method that check is a point is in the MSET (you probably want to skip the last bruteforce check  ).
It can be optimized by using some temporary variable if you're computing on a CPU and not a GPU. (optimization that doesn't work on GPU because of read-after-write register latency of ~24 cycle)

Code:

bool isInMSet(
    float cr,
    float ci,
    const uint maxIter,
    const float escapeOrbit)
{
    int iter = 0;
    float zr = 0.0;
    float zi = 0.0;
    float ci2 = ci*ci;
    float temp;

    //Quick rejection check if c is in 2nd order period bulb
    if( (cr+1.0) * (cr+1.0) + ci2 < 0.0625) return true;

    //Quick rejection check if c is in main cardioid
    float q = (cr-0.25)*(cr-0.25) + ci2;
    if( q*(q+(cr-0.25)) < 0.25*ci2) return true; 


    // test for the smaller bulb left of the period-2 bulb
    if (( ((cr+1.309)*(cr+1.309)) + ci*ci) < 0.00345) return true;

    // check for the smaller bulbs on top and bottom of the cardioid
    if ((((cr+0.125)*(cr+0.125)) + (ci-0.744)*(ci-0.744)) < 0.0088) return true;
    if ((((cr+0.125)*(cr+0.125)) + (ci+0.744)*(ci+0.744)) < 0.0088) return true;

     //BRUTEFORCE CHECK
    while( (iter < maxIter) && ((zr*zr+zi*zi) < escapeOrbit) )
    {
        temp = zr * zi;
        zr = zr*zr - zi*zi + cr;
        zi = temp + temp + ci;
        iter++;
    }

    if( iter < maxIter)
    {
        return false;
    } else {
        return true;
    }
    //END OF BRUTEFORCE CHECK

}

[ker2x]

and alternative openCL implementation :

Code:

//Check if choosen point is in MSet
bool isInMSet(
    const float2 c,
    const uint minIter,
    const uint maxIter,
    const float escapeOrbit)
{
    int iter = 0;
    float2 z = 0.0;

    if( !(((c.x-0.25)*(c.x-0.25) + (c.y * c.y))*(((c.x-0.25)*(c.x-0.25) + (c.y * c.y))+(c.x-0.25)) < 0.25* c.y * c.y))  //main cardioid
    {
        if( !((c.x+1.0) * (c.x+1.0) + (c.y * c.y) < 0.0625))            //2nd order period bulb
        {
            if (!(( ((c.x+1.309)*(c.x+1.309)) + c.y*c.y) < 0.00345))    //smaller bulb left of the period-2 bulb
            {
                if (!((((c.x+0.125)*(c.x+0.125)) + (c.y-0.744)*(c.y-0.744)) < 0.0088))      // smaller bulb bottom of the main cardioid
                {
                    if (!((((c.x+0.125)*(c.x+0.125)) + (c.y+0.744)*(c.y+0.744)) < 0.0088))  //smaller bulb top of the main cardioid
                    {
                        while( (iter < maxIter) && (z.x*z.x + z.y*z.y < escapeOrbit) )      //Bruteforce check  
                        {
                            z = (float2)(z.x * z.x - z.y * z.y, (z.x * z.y * 2.0)) + c;
                            iter++;
                        }
                        if( (iter > minIter) && (iter < maxIter))
                        {
                            return false;
                        }
                    }
                }
            }
        }
    }
    return true;
}

[fractalmind]

In my opinion you only know which mandelbrot algorithm is really optimized (that is in real life) by benchmarking some deep zooms. Otherwise this discussion is purely academical.

You can only use double precision (64 bit)  accuracy with GPU unless you write some arbitrary precision library (like GMP) in opencl (if that is possible) or similar gpu api.

[ker2x]

In my opinion you only know which mandelbrot algorithm is really optimized (that is in real life) by benchmarking some deep zooms. Otherwise this discussion is purely academical.

But it's very usefull for the buddhabrot (my opencl code above was for my older buddhabrot renderer).

[fractalmind]

Well buddhabrot is a little different creature...

I have been doing some benchmarking against Ultra Fractal with my program.

For example:

deep2.fdz (http://personal.inet.fi/koti/fractalus/fractalus570.zip in gallery folder)

center coordinates:-1.744416922012436877855459748238956035814664967683655064429/2.2045963513788370573266571752910479871072765644643956489469E-2 magnification:8.4504202435435706944886207459268813904116214019983922892912E44

Max iterations: 10 000 bailout=512 precision: 176 bits (50 decimals). Image size 640x480

Ultra fractal 3min 20 s(with guessing, periodicity checking normal) using 8 threads

Fractalus (my program) 2 min 52s also using 8 threads with 32bit gmp.dll

Machine: 2600k oveclocked to 4600mhz

I use slightly modified algorithm from my previous post in this thread. No special optimizations (guessing or "cheating" since they do not seem to make much difference especially with deeper magnifications.)

[PurpleBlu3s]

When I made the thread I was actually thinking about the Buddhabrot, and that's what I'm using it for.

[ker2x]

When I made the thread I was actually thinking about the Buddhabrot, and that's what I'm using it for.

  

[Duncan C]

In my opinion, cardioid checking (and testing for the low order buds) isn't a very useful strategy for optimising the calculation of the Mandelbrot fractal. This technique is only good for very low zooms, which only require a low number of iterations and are therefore fast without any optimisation.

More effective techniques for optimising the calculation (that work at higher zoom levels) are periodicity checking (also known as obit detection) and adjacency optimisation. I recommend having a read of Robert Munafo's Encyclopedia of the Mandelbrot Set: http://mrob.com/pub/muency.html . He describes various techniques for adjacency optimisation: Solid guessing (successive refinement), boundary tracing and Mariani / Silver (also known as tesseral). All these techniques will help speed up the calculation inside the set (the black area).

There are two approaches to periodicity checking, one of which is described on Robert's website (the RHOP algorithm, the other being the binary back off algorithm). Both techniques will add overhead to the calculation but will greatly speed up detecting points inside the set. I believe the RHOP algorithm is faster than the BBO algorithm, but adds greater overhead (could be wrong though). I use the latter algorithm in my Mandelbrot program.

I have also found this page with a listing in Java, which you may find useful (includes cardioid / bulb detection and period checking - using the binary back off algorithm):

http://en.wikipedia.org/wiki/User:Simpsons_contributor/periodicity_checking

zenzero-2001

I agree with zenzero about the cardioid optimization being of very limited usefulness. I did a lot of research on various heuristics for Mandelbrot sets, and finally settled on using a boundary following algorithm.

The Mandelbrot set is connected, and so are the iteration bands around it. The connections are infinitely thin, and not detectable at any given scale, but you can be sure that "blobs" of Mandelbrot points will never contain non-mandelbrot points inside. You can therefore walk the outside of any group of Mandelbrot pictures, always turning in the same direction (a common maze solving algorithm) Then, once you get back to your starting point, you can flood-fill the inside.

This also works for any iteration band surrounding the Mandelbrot set. The deeper the zoom the less useful boundary following tends to be for non-Mandelbrot points, but it is even more useful for Mandelbrot points because you typically need to increase the max iteration count for deeper zooms, so points that are members of the Mandelbrot set become very expensive to calculate.

Boundary following works for ANY portion of the Mandelbrot set (with varying degrees of success), whereas plotting known shapes like the main cardioid , "Baby Mandelbrot" cardioid", and main disks only works for zooms that include those small numbers of shapes.

My application, FractalWorks, includes boundary following as well as multi-threading and symmetry calculations. I even got boundary following working for Julia sets. I found through testing that even for "Fatou Dust" plots (disconnected Julia set plots who's origin point is not a Mandelbrot point) I can create perfect plots without any holes if I start boundary following above 10 or 15 iterations.

I got the idea for boundary following from Dr. Evgeny Demidov, who wrote a VERY impressive Mandelbrot and Julia set renderers for the web in Java that was dramatically faster than early versions of FractalWorks. He published the source code for his Java program, which gave me enough information to implement my own boundary following algorithm.
FractalWorks has preferences that let you turn boundary following on and off, and tracks the total pixels and iterations calculated vs. those skipped because of boundary following. For plots that have large blobs of Mandelbrot points, the effect can be quite dramatic. It can give 10-fold or better improvements in speed.

If you request Distance Estimate or fractional iteration values for a plot, the program can't use boundary following for any points except Mandelbrot/Julia set points, but that is still a big win for any plot that has any "Baby Mandelbrots" in it.

FractalWorks doesn't support orbit trap coloring of Mandelbrot points, so I can always use boundary following for those points. If you wanted to use orbit trap calculations you would probably have to brute-force Mandelbrot points.

FractalWorks is only available for Macintoshes. It's available in the Mac App Store (MAS) if you're interested in trying it out. Here's the link:

FractalWorks in the Mac App Store

FractalWorks requires Mac OS 10.6.6 or later.
-----
EDIT:
I did a benchmark plot of the whole Mandelbrot set (real width 3.0) at 10,000 max iterations.

With boundary following and symmetry turned off, it took 7.583 seconds, and required calculating 1,089,686,150 iterations.

With boundary following turned on (and symmetry still turned off) the same plot takes between 0.25 seconds and and 0.367 seconds. (At plot times that short, things like screen refreshing and when the app gets swapped into the background have a disproportionately large effect on plot times, so the recorded plot times can vary by a fair margin.). The plot  required calculating 59,855,336 iterations, or about 5.5% of the total iterations for the brute force method.

Zooms into the cleft of Seahorse Valley or one of the cardioids can have even more dramatic speedups. I just did a plot of the cleft of the main cardioid that's mostly black with boundary following turned on, and FractalWorks is able to skip 98.93% of the total iterations that would be required for the "brute force" approach. If I were to turn symmetry calculations back on for that plot the speed roughly doubles again. These speedups work for the main body of the Mandelbrot set, or for any of the "Baby Mandelbrots"