Understanding Perturbation

[nitroxis]

Hello.

I've been working on a simple Mandelbrot zoomer using OpenGL Compute Shaders (a relatively new kind of GPU programs) and multi-precision fixed-point numbers. The program itself works pretty well and is faster than most CPU implementations. But as with all implementations of the standard Z=Z*Z+C technique, the performance depends heavily on the magnification (i.e. precision of the fixed-point numbers) and number of iterations used.

I recently came across the perturbation technique implemented in SuperFractalThing and others and need some help understanding and implementing it. I tried to understand the source code of SuperFractalThing, but it's rather hard for me to get an idea of what's going on.

So the basic idea of using perturbation is that the image (except for the reference points) can be calculated using hardware floating point numbers (which sounds great, considering that's what GPUs are best at). The first thing to do is choosing a reference point and iterating it using arbitrary precision numbers. Then the rest of the image can be calculated using delta iteration based on that reference orbit. That's what I don't really get.

First, the reference orbit is in arbitrary precision. Do I just convert it to floating-point numbers in order to use it in the delta iteration? Second, depending on how far you've zoomed in, the delta (i.e. the distance between two pixels) get extremely small, so small that they can't be reasonably represented using floating-point numbers. How does that work? Can I somehow scale up the deltas in order to work with the larger numbers? Third, how do I find out when the iterated value's magnitude is larger than the bailout value? I only have the deltas, so I suppose I have to add it to the reference point and check the magnitude of reference[n]+delta[n] > bailout, but then again, the delta is supposed to be a floating point number whereas the reference orbit is arbitrary precision and not necessarily small, so a lot of precision would be lost when adding the two.

I'd be glad if someone familiar with the technique could explain a little more how exactly the "low precision" part is supposed to work.

[youhn]

Use normal floating points down to 1x10−308
Use double precision for everything down to 3.65×10−4951

Bailout condition need not be so high. Kalles Fraktaler uses "only" 100 for smooth method, 2 for normal.

I cannot help you further, since I do not understand the computational steps in detail for the perturbation thing.

[nitroxis]

I'm using the "standard" bailout of |z| > 2. I was talking about adding the reference point (usually a number between -2 and 2) to the delta (a number that is extremely small for high levels of magnification) and the loss of precision in this step.

[Kalles Fraktaler]

First, the reference orbit is in arbitrary precision. Do I just convert it to floating-point numbers in order to use it in the delta iteration?

Yes!

Second, depending on how far you've zoomed in, the delta (i.e. the distance between two pixels) get extremely small, so small that they can't be reasonably represented using floating-point numbers. How does that work?

You can with the double datatype represent values as small as 1e-308.

Can I somehow scale up the deltas in order to work with the larger numbers?

You can multiply the high precision values before converting them to double values with up to 1e308. You then need to divide the values with the same scaling value.
By doing so, you can use ordinary double datatypes up to 1e616, in theory.

In practice, Kalles Fraktaler is using scaling from 1e301 to 1e600 and then the long double datatype is used from 1e600 to 1e4900. And beyond that I use a custom datatype with double as mantissa and integer as exponent. See floatexp in the source of Kalles Fraktaler. (but you could of course use scaled long double to 1e9800)

Third, how do I find out when the iterated value's magnitude is larger than the bailout value? I only have the deltas, so I suppose I have to add it to the reference point and check the magnitude of reference[n]+delta[n] > bailout, but then again, the delta is supposed to be a floating point number whereas the reference orbit is arbitrary precision and not necessarily small, so a lot of precision would be lost when adding the two.

You have already converted the arbitrary precision values to doubles, so reference[n]+delta[n] grows beyond the bailout value. Also the delta value grows from extremely small to a value comparable to the bailout value in magnitude.

And then you will find out that the low precision prevents some areas of a view to be rendered correctly and that glitches will occur in your images. Some of these glitches are one-colored blobs, some with incorrect structures. Then you are ready for the next step, to use Pauldelbrot's glitch detection method to detect these pixels, which need additional high precision references to be rendered correctly.

But you will have a lot of fun!

And don't forget Series Approximation, that allows you to skip a lot of iterations and that is many times even more time saving than what you get from using hardware datatypes instead of high precision.

[nitroxis]

Thank you, I've got a better idea of how it works now.

Yes, I already know about glitches and the like. It was more about the fundamental ideas of the algorithm that I didn't fully understand.
I got a somewhat working implementation now (which just uses the center of the screen as reference point).

[cKleinhuis]

kalles pure cpu implementation already is awesomly fast, i am looking forward to your attempt using a gpu based implementation, its kinda interesting how it develops, some people that implemented it said they use more than one gigabyte to store such reference high precision variable, it comes down to pure memory to get deeper, i have 16 gigs available here, lol i would love uising 8 or more gigs to get a view of regions of the mandelbrot where no one has gone before

[nitroxis]

Well, it would require the RAM of the graphics card, not the system RAM. But gigabytes sound a bit much, considering you only need 16 bytes for one complex number (2 doubles). It would mainly depend on the number of iterations used. One GB of VRAM would already be enough to hold 67108864 complex numbers. That'd be enough for millions of iterations, depending on whether you use series approximation, the derivatives and so on. The memory required for the output would only depend on the image size used to render the fractal. So I guess it would scale arguably well, but I'll have to see about that.

One thing to reduce the memory footprint of the application would be to only calculate say 1000 iterations of the reference point, then calculate pixels around it with delta iteration, calculate the next 1000 iterations of the reference point, continue with the delta iterations and so on. This way one could achieve good CPU and GPU parallelism, as the GPU can compute the pixels around the reference point while the point itself is still being calculated by the CPU.

[cKleinhuis]

right, but we are talking about depths of a billion or even more iterations!

[3dickulus]

re:storage this image, "Sircle" by dinkydau, http://www.fractalforums.com/index.php?topic=18254.msg76549#msg76549 took over 300M to store the reference data at almost 13,000,000 iterations. all of the stored data is double (64bit) format. by this one could infer that 1G will hold about 40,000,000+ iterations worth of data.  well, at least in the way SFTC stores it, mX[], mXi[], mDistance_to_edge_sqr[]. With series approx I think this could be cut down considerably by not storing the "skipped" iterations and just maintaining the current working range required for the image.

edit:"Sircle" location from dinkydau, image rendered by SFTC

[Alef]

What is deepest image so far? Hadn't followed this stuff for a while? That example is kind of nice throught it isn't 800x600 at least.

[Kalles Fraktaler]

Here is the deepest image(s)

<a href="http://www.youtube.com/v/WgzW_jr0xSw&rel=1&fs=1&hd=1" target="_blank">http://www.youtube.com/v/WgzW_jr0xSw&rel=1&fs=1&hd=1</a>

[cKleinhuis]

Here is the deepest image(s)

<a href="http://www.youtube.com/v/WgzW_jr0xSw&rel=1&fs=1&hd=1" target="_blank">http://www.youtube.com/v/WgzW_jr0xSw&rel=1&fs=1&hd=1</a>

lol, its just a video
i know you have been deeper, just did not save and publish

[Kalles Fraktaler]

lol, its just a video
i know you have been deeper, just did not save and publish

No, I have never been deeper!
The "Find Minibrot" function was fortunately successful on this location and was working for about a month!

If you use KF with standard resolution 640x360, zoom lever 32 is about the largest to have some control of where to zoom.
Zooming manually to 1e10000 would require 6644 clicks!
I'll probably never go deeper than this...

[nitroxis]

Alright, I got my program working with plain delta iteration. It's able to recognize glitches, however it's not correcting them yet. It currently uses the center of the screen as a reference point. This is not optimal, but I don't really know what to do here. How do I find a relatively good reference point to start with (given the position in the complex plane and a zoom factor)? After the first pass is rendered, do I just use the center of the glitch blobs as new reference point or is there a better way to do this? I don't really understand how SFT is doing it.

Here's what I got so far (red means glitch).

[Kalles Fraktaler]

Is K.I.Martin's SFT solving glitches at all?