Perturbation formula for burning ship (hopefully correct :P)

[laser blaster]

I took a shot at deriving the perturbation formula for the burning ship fractal. I need someone to check it for correctness though, because I'm not really a math guru.

So with the burning ship formula, the component-wise abs function means you need to split it into real and imaginary components. Let Z be the reference point, and let r and i be the real and imaginary components of Z.

Z_n+1:    
    r_n+1 = r_n² - i_n² + r₀
    i_n+1 = 2 * abs(r_n * i_n) + i₀
(I simplified out some of the abs()s, the result will be the same)

Let W be the pixel we are going to iterate, and let x and y be the real and imaginary components of W.
Let D be the difference between the pixel point and reference point, let a and b be the real and imag parts of D

W₀ = Z₀ + D₀  
D₀ = W₀ - Z₀

D_n+1 = W_n+1 - Z_n+1

D_n+1:
    a_n+1 = (x_n² - y_n² + x₀) - (r_n² -i_n² + r₀)
   edit:
    b_n+1 = 2*(abs(x_ny_n) + y₀ - abs(r_ni_n)-i₀)
   b_n+1 = 2*abs(x_ny_n) + y₀ - 2*abs(r_ni_n)-i₀

Now, putting D_n+1 in terms of r_n,i_n,a_n, and b_n only (getting rid of all the x's and y's):

a_n+1 = ( (r_n + a_n)² -  (i_n + b_n)² + r₀ + a₀) ) - (r_n²  - i_n² + r₀)
    edit:
    = (r_n² + 2a_nr_n + a_n² - i_n² + 2i_nb_n + b_n² +r₀ + a₀) - (r_n² - i_n² + r₀)
    = (r_n² + 2a_nr_n + a_n² - i_n² - 2i_nb_n - b_n² +r₀ + a₀) - (r_n² - i_n² + r₀)
    some of the terms cancel...
    = 2a_nr_n + a_n² - 2i_nb_n - b_n² + a₀
    Now all the terms in the equation for a are "small" when D is small.

b is harder to simplify due to the non-removable abs() function. But here goes: (continued next post)

[laser blaster]

edit:
b_n+1 = 2( abs((r_n+a_n)*(i_n+b_n)) + i₀ + b₀ - abs(r_ni_n) - i₀ )
    = 2( abs(r_ni_n+r_nb_n+a_ni_n+a_nb_n) + b₀ - abs(r_ni_n) )

b_n+1 = 2*abs((r_n+a_n)*(i_n+b_n)) + i₀ + b₀ - 2*abs(r_ni_n) - i₀
    = 2( abs(r_ni_n+r_nb_n+a_ni_n+a_nb_n) - abs(r_ni_n) ) + b₀


Now the important thing to note about the above formula is that 3 of the terms inside the first abs() contain a_n or b_n as a factor, and so those three terms will be small when D is small. There is one term, r_ni_n, that doesn't have a or b as a factor (so it's too "large"), but there's a tantalizing  "- abs(r_ni_n)" after that that looks like it could be made to cancel with the earlier term. However, the  first occurence of r_ni_n is wrapped up in an abs() function with other terms, which complicates things.

To understand how to make the large terms vanish, I simplified the problem to it's basics: given the formula abs(c+d) - abs(c), try to find the answer without having the result depend on c. (c stands in for r_ni_n (the "large" term), and d stands in for the sum of the other, smaller terms) (comparisons with c in conditional statements are fine, as you can use a rounded double precision representation of c and experience no precision loss, as long as your other comparison operand is "small" enough to be represented as a double(d meets this criteria)).


The result can be expressed without an abs, using conditionals:

abs(c+d) - abs(c):
    if c>0:
        if c+d > 0 (e.g. d > -c):
            result = d
        else if d == -c:
            result = d
        else if d < -c:
            result = -d -2c
    else if c==0:
        result = abs(d)
    else if c < 0:
        if c+d>0 (e.g. d > -c)
            result = d + 2c
        else if d == -c:
            result = -d
        else if d < -c:
            result = -d

That should give the proper result for all c and d, I think (I hope).
Notice that I was able to express the result in terms of d alone (except for the comparisons between d and c, but as I mentioned above, these can be done exactly in the precision of the least precise operand, so we can use fast native doubles), except for two cases, where the result still involves a (-d - 2c, and d + 2c). Upon first glance it seems that we’re screwed, because c is too “large” for the perturbation technique to work. But notice that the result only involves c in when either: c is positive and d is less than -c, or c is negative and d is greater than -c. So in either case, the magnitude of d  (the supposed “small” term) is already LARGER than the magnitude of c, and what’s more, the two have opposite signs. So the addition of the 2c or -2c term will not result in any loss of significant digits of the difference between adjacent pixel points. And, as I learned in Pauldelbrot’s glitch-fixing thread, the magnitude of the D (and here I mean delta, the difference between the reference and pixel points) term isn’t what’s important, the important thing is that there be enough bits of precision to accurately represent the difference between adjacent points, and as long as that’s the case, there won’t be any glitches.
    
Now just replace d with r_nb_n+a_ni_n+a_nb_n and c with r_ni_n, calculate the result using the above algorithm with the conditionals, multiply that by 2, and add to that b₀,  and you’ve got b_n+1. Now you can find the real and imaginary components of D_n+1.

The one thing I didn't attempt was series approximation. I have a feeling it's possible, but it will obviously be two real-valued functions, each in two varaibles, plus with special considerations for the abs() function, rather than complex-valued polynomial function of a complex variable. Much more complicated than the Mandy approximation

I might have made a typo just typing in my work on paper, so I will check for that tomorrow.

[laser blaster]

Uh oh, I think I see a mistake.  

I need to learn to check before posting.  Looks like I have to redo the work for b. I think it should still be possible.

[DarkBeam]

Uh-oh, abs() makes bad jokes when you derive... Good luck

[Kalles Fraktaler]

Cool, I hope you succeed finding the formula!

I guess Series Approximation is not possible since it is based on derivatives...?

[laser blaster]

Okay, the abs() was throwing me for a loop... but it is valid to replace abs(x)*abs(y) with abs(x*y), then substitute x any y with an equivalent expression, and expand out the multiplication, right? So that part was fine. But I made another mistake, which I'm gonna fix now.

Okay, I fixed a few things. It might be correct now, but I'll check over it a few more times to make sure it all works out.

[laser blaster]

Wooo! It works!

I wrote a quick test in Fragmentarium. I calculated the reference point in double precision (of course, due to the nature of pixel shaders, I have to recalculate the reference point for every pixel, so it runs slower than just using double precision in the first place, but a CPU-based implementation will have no such limitations), and calculated delta in just single precision, and now I can zoom to double precision levels without pixelization.

I attached a comparison picture showing the difference without perturbation (It'll be up in a sec)

I also attached the .frag file along with a slightly changed version of the Progressive2D.frag include file ( I just added 1 line of code to enable native doubles). Note that I actually calculate the reference point in the main loop along with the delta. Trust me, it would work just as well if the reference point were precomputed and stored in a array.

[laser blaster]

Here's the comparison pic.

[Kalles Fraktaler]

Wooohooo! It works!

But your formula

b_n+1 = 2( abs(r_ni_n+r_nb_n+a_ni_n+a_nb_n) - abs(r_ni_n) ) + b₀

don't work since r_nb_n+a_ni_n+a_nb_n can be much smaller/bigger than r_ni_n when zooming deeper than e18.
So I had to use your conditions...

[laser blaster]

Awesome! I can't wait to try it!

[Kalles Fraktaler]

I wish I could use the "Power" field to set an arbitrary power also for the Burning Ship fractal in Kalles Fraktaler.
That would be so cool!

However, the absolute values complicates things a lot.
I made a simple complex class for the higher powers of Mandelbrot, that only implemented addition, subtraction and multiplication. I use it from power 4, because it get too complicated to expand all the terms and too easy to make mistakes.

When making full precision calculations of a pixel in the Burning Ship fractal, I have found out that the standard Mandelbrot formula can be used with the only addition to set the real and imaginary parts to absolute values for each iteration (and the square (2) can be replaced by any value):
while(!bailout){
  X = X² + X₀
  X.r = abs(X.r)
  X.i = abs(X.i)
}


Unfortunately this does not work on the perturbation formula that calculates from a reference:
while(!bailout){
  D = 2*X*D + D² + D₀
  D.r = abs(D.r)
  D.i = abs(D.i)
}

Firstly, it seems the absolutes needs to be applied considering how the real and imaginary parts are grouped in the expanded expressions.
Secondly, there are these "abs(c+d) - abs(c)" situations, which in higher depths (>e18) don't give accurate result and must be solved by the conditions that laser blaster created.

Maybe there is a way to generalize arbitrary powers of Burning Ship, but it seems very complicated. I did an attempt to the 3rd power (using the same variable names as laser blaster), but I haven't tried to program it yet:

r = abs(r³) - abs(3*r*i²) + r₀
i = abs(3*r²*i) - abs(i³) + i₀

a = abs(x³) - abs(3*x*y²) + x₀ - (abs(r³) - abs(3*r*i²) + r₀)
b = abs(3*x²*y) - abs(y³) + y₀ - (abs(3*r²*i) - abs(i³) + i₀)

a = abs((r+a)³) - abs(3*(r+a)*(i+b)²) + r₀ + a₀ - abs(r³) + abs(3*r*i²) - r₀
a = abs((r³ + 3*r²*a + 3*r*a²+a³)) - abs(3*(r+a)*(i²+2*i*b+b²)) + a₀ - abs(r³) + abs(3*r*i²)
a = abs((r³ + 3*r²*a + 3*r*a²+a³)) - abs(r³) + abs(3*r*i²) - abs(3*r*i² + 6*r*i*b + 3*r*b² + 3*a*i² + 6*a*i*b + 3*a*b²) + a₀
...we have two abs(c+d) - abs(c) situations here (one is reversed)

b = abs(3*(r+a)²*(i+b)) - abs((i+b)³) + (i₀+b₀) - abs(3*r²*i) + abs(i³) - i₀
b = abs(3*r²*i+6*r*a*i+3*a²*i + 3*r²*b+6*r*a*b+3*a²*b) - abs(i³+3*i²*b+3*i*b²+b³) + b₀ - abs(3*r²*i) + abs(i³)
b = abs(3*r²*i+6*r*a*i+3*a²*i + 3*r²*b+6*r*a*b+3*a²*b) - abs(3*r²*i) + abs(i³) - abs(i³+3*i²*b+3*i*b²+b³) + b₀
...and we also have two abs(c+d) - abs(c) situations here, again one reversed

So I guess it is necessary to expand everything for each power and identify all abs(c+d) - abs(c) situations and grouping, and the 4th power would be enormously complicated...
Unless the complex class could be modified to be able to handle the absolutes automatically.
But I don't know how...

[Kalles Fraktaler]

Unfortunately my attempt was not successful... 

[stardust4ever]

I'm bumping this thread because I would like to see it done. The third order Burning Ship is one of the most beautifully symmetrical fractals I have seen, with three needles. It was originally popularized by a video posted by HPDZ.

The code I developed in 2012 for the floating point portion of the FX plugin works like this:

Code:

Cubic Burning Ship:

		zi = ((zrsqr * 3.0) - zisqr) * abs(zi) + JuliaI;
		zr = (zrsqr - (zisqr * 3.0)) * abs(zr) + JuliaR;
		zisqr = zi * zi;
		zrsqr = zr * zr;

Cubic Buffalo:

		zi = abs(((zrsqr * 3.0) - zisqr) * zi) + JuliaI;
		zr = abs((zrsqr - (zisqr * 3.0)) * zr) + JuliaR;
		zisqr = zi * zi;
		zrsqr = zr * zr;

Although we had several missed attempts at getting the 3rd order formula right, compare to these 2nd order formula that Panzerboy initially created:

Code:

Burning ship:

		zi = abs(zr * zi) * -2.0 + JuliaI;
		zr = zrsqr - zisqr + JuliaR;
		zisqr = zi * zi;
		zrsqr = zr * zr;

Celtic:

		zi = zr * zi * 2.0 + JuliaI;
		zr = abs(zrsqr - zisqr) + JuliaR;
		zisqr = zi * zi;
		zrsqr = zr * zr;

Buffalo:

		zi = abs(zr * zi) * -2.0 + JuliaI;
		zr = abs(zrsqr - zisqr) + JuliaR;
		zisqr = zi * zi;
		zrsqr = zr * zr;

As seen above, the 2nd order abs fractals can be generated by applying the abs() command to the entire real or imaginary part of the equation.

In fact, these are all subsets of a larger group of twelve 2nd order Mandelbrot variant fractals I have categorized, including the Mandelbar and Perpendicular sets. Ten of those formulas including the standard M-set are made available in Fractal Extreme as a set of plugins by Panzerboy, some of which I helped create. The two fractal formulas that got omitted are nothing to write home about. I'll save further discussion regarding the twelve Mandelbrot abs() variants, as well as some of my recent musings into the 3rd order, for another topic.

Applying the abs() command only to the entire Zr or Zi after exponent does not get the same results for higher orders. In second order, if both the Zi and Zr components are absoluted prior to squaring, we get essentially the same result as applying the abs() part on the imaginary side post squaring. The reason for this is that the real component contains only square parts, making an abs() operation prior to squaring moot. The Imaginary component contains the product of Zr and Zi, so the result will be positive whether the abs() command is applied before or after. The "-2" multiplicand is only to flip the imaginary channel, making the ship appear mast up as it is most commonly rendered. For the third order version of the Burning Ship, there is only one operand per coordinate that needs the abs() function, due to distributive property of mathematics.

I really don't have a grasp of the actual mathematics behind perturbation theory, but if it works with the burning ship fractal, it should work with the 3rd order as well? It seems that the pixels mostly stick together as their values orbit the set, and the change (delta) values between adjacent pixels changes little during the course of the journey, enabling lower precision than the main orbit. Do the abs() commands really screw it up that bad? It seems the base fractal would not be possible if it did.

Honestly I would be happy with just the 2nd and 3rd order version. I have viewed the 4th and higher higher order versions of the Burning Ship fractals and on first glance (using slow as molassis UF5), they simply don't have the ornate symmetry of the 2nd and 3rd order versions.

[youhn]

That abs operation makes a very big difference between the mandelbrot set and the burning ship. This addition makes the burning ship non-analytical. If I got it right, this should mean that the perturbation thing should NOT work on the bs fractal. But somehow it does. Perhaps the reason is what you just wrote:

The reason for this is that the real component contains only square parts, making an abs() operation prior to squaring moot. The Imaginary component contains the product of Zr and Zi, so the result will be positive whether the abs() command is applied before or after.

If this is the case, the perturbation thing should work on anything with even powers.

[stardust4ever]

That abs operation makes a very big difference between the mandelbrot set and the burning ship. This addition makes the burning ship non-analytical. If I got it right, this should mean that the perturbation thing should NOT work on the bs fractal. But somehow it does. Perhaps the reason is what you just wrote:

If this is the case, the perturbation thing should work on anything with even powers.

So abs() operands are allowed within the imaginary part but not the real part? Or is it a requirement the abs operation needs to be performed after the exponent (the entire Zi or Zr, not just a portion of it)? One way to look at it is like a fast driving car. Mandelbrot is fully connected like a highway with gradual, wide turns. Many shapes in the abs sets are extremely angular and disconnected. No matter how well you think your car handles, in real life a Dukes of Hazzard style jump or a hairpin turn going 155mph will cause you to crash... If a pixel encounters an abs operation prior to the bailout (for instance the delta orbit of the pixel gets a sign change at some point compared to the reference orbit resulting in a detour of the orbit path, which is probably likely at some point given the frequencies of angular and reflected symmetries within the fractal), does this automatically cause the float point estimator to fail?

If so, then why does the quadratic BS work?