Welcome to Fractal Forums

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. :D

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)

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.

Uh oh, I think I see a mistake. :o  :sad1: :'(

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.

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

Cool, I hope you succeed finding the formula!

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

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.

Wooo! It works! :D :joy:

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.

Here's the comparison pic.

Wooohooo! It works! :D :joy:

But your formula
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... :)

Awesome! I can't wait to try it! O0

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...

Unfortunately my attempt was not successful...  :hurt:

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:


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

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.

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? :nerd:

Because of laser blaster's excellent if-conditions, which I don't even fully understand

(because my expertise is not to create the fastest most optimized or the most beautiful easy-to-understand code. My expertise is to very quickly implement an idea, often someone else's, into a working program and that is also what I do for a living :) )

Perturbation theory isn't limited to analytic functions.It's much more convenient to work out for analytic functions, as you can use complex arithmetic for those, and you can operate on complex numbers the same way as you do for reals. However, when you have a non-analytic function, the complex numbers no longer apply, and you have to treat it as a function of two real variables, which means more work.

The abs() folding makes it even trickier to get accurate results. This is because you can't easily make the large terms cancel: How can you simplify abs(a+c) - abs(a) to an expression involving c only? Well, you can't. However, under certain conditions, the function can be expressed entirely without abs. For example, if c and a are both greater than 0, then abs(a+c) = a+c and abs(a) = a. In this way, by analyzing the function case-by-case, you can simplify and make the large terms cancel for most cases. However, there are still some cases where you can't cancel the large terms. But it just so happens that in all such cases, the delta term (the perturbed orbit) is larger than any of the terms from the reference orbit(which are usually larger),  and the two terms have opposite signs, so adding in the typically "large" term from the reference orbit actually brings the sum closer to zero, so no precision is lost.

I'm looking into deriving the cubic perturbation formula. I've made some progress, but I'm not yet sure if it'll all work out. On the other hand, for powers that are powers of two (2,4,8,etc), deriving the perturbation formula should be trivial.

Sweet. So you have to apply a bunch of conditional logic (If Then, etc...) to a variety of formula variants in order make each abs() function work? Sounds complicated.

For Quadratic Burning Ship,

*Burning ship:


You can apply the abs() function to Zi0 and Zr0 individually prior to iteration, or apply the abs() function once to (zr*zi) because the result will be identical either way.

For Cubic Burning Ship,


applying the abs() to just zi1 will not work and result in an erroneous fractal. Applying the abs() to zr0 and zi0 respectively prior to the equation is the same net effect as applying the abs() to zi and zr sides of the equation individually. This is again due to the fact that X^2 is always positive.

If applying the perturbation formula to the cubic fractal requires two abs() commands instead of one as with the quadratic, then there are four conditional equations to be selected based on signage, as opposed to two with the burning ship. I guess that does complicate things a bit, but it should make the fractal impossible to calculate.

Thanks for the effort in working on this! O0

That sounds really exiting. I keep my fingers crossed :pray2:

With stardust4ever's formula, together with laser blaster's if-conditions, the cubic burning ship is actually as "easy" as the ordinary burning ship.

Sometimes, it is enough encouraging just to read that someone writes it is possible. So when laser blaster says it's possible, I believed it.

But I don't think I have the time to make a it available in KF for you all until next week, unfortunately...

(http://i717.photobucket.com/albums/ww173/prestonjjrtr/Smileys/worship.gif)

 :cantor_dance:

 O0 O0 O0

https://www.youtube.com/watch?v=_rL6pI8-ftU

No, thank you! Awesome, awesome video btw...
 :thanks2: :toast: :beer: :beer: :beer:

Thanks

Don't forget to check out Kalles Fraktaler 2.7 :)

BTW stardust, do you have a formula for cubic celtic, in the same easily accessible way as you have provided the others here in this thread?
(http://www.fractalforums.com/index.php?action=dlattach;topic=2481.0;attach=6861;image)
I had a flow yesterday evening and implemented 4 formulas without a single failure. Only the cubic celtic is missing ;)

Expect a PM shortly.

I am working on a database of fractal formulas of base 2 and 3 that add abs() formulas to the fractal, as well as negating the imaginary side. I took a break for a bit due to recent surgery but I am recovered and at home now. I have cataloged a total of 12 fractals in second power and am now working on 3rd power.

In the second power fractal, I found 8 unique symmetrical fractals and 8 chiral pairs of asymmetric fractals. Because I don't count mirrors or rotations, that makes 12 unique second order formula. These fractals were fabricated by applying the abs() function to certain key components of the equation, or not, and optionally sign-changing the imaginary part of the equation. Applying a sign change to the real portion always produces a duplicate of one of the other fractals, with the needle pointing East instead of West, so that eliminates most of the duplicates.

The way I plan on slitting the 3rd order equation, will result in 64 distinct combinations, however many of them will also be rotations or mirrors.

I will eventually publish both lists, but for the time being, I can send you a copy of the 2nd order list (with footnotes) and some of the entries discovered so far in the 3rd.

As for deep fractal zooming of the cubic variety, the Cubic Burning Ship and the Cubic Buffalo (plus a new 3rd order hybrid Burning Ship/ Buffalo I discovered which takes the imaginary portion of one and the real portion of the other) seem to be the most promising.

I found out that the Cubic Buffalo may not be possible to render with perturbation on certain circumstances.
Locations very close to - but not on - the x-axis are not rendered correctly. Different references create different views - however good looking and not the ordinary pixelated glitches that normally occur due to too low precision. Instead, the detected glitches are filled with incorrect pattern, and if I put the main reference somewhere else the view is not the same.
Any help is much appreciated, even though I understand it is very difficult...
Example:

Here is my calculations, from stardust4ever's easily accessible formula:
x-axis:
zr = abs((zrsqr - (zisqr * 3.0)) * zr) + JuliaR;
abs((x*x - (y*y * 3.0)) * x) + x0  -  abs((r*r - (i*i * 3.0)) * r) - r0;
abs(((r+a)*(r+a) - (y*y * 3.0)) * (r+a)) + a0 - abs((r*r - (i*i * 3.0)) * r);
abs( (r*r + 2*r*a + a*a - 3*i*i - 6*i*b - 3*b*b) * (r+a) ) - abs(r*r*r - 3*r*i*i) + a0;
abs( r*r*r-3*r*i*i + 2*r*r*a + r*a*a - 6*r*i*b - 3*r*b*b + r*r*a + 2*r*a*a + a*a*a - 3*i*i*a - 6*i*b*a - 3*b*b*a ) - abs(r*r*r-3*r*i*i) + a0;
abs( r*r*r-3*r*i*i + 3*r*r*a + 3*r*a*a - 6*r*i*b - 3*r*b*b + a*a*a - 3*i*i*a - 6*i*b*a - 3*b*b*a ) - abs(r*r*r-3*r*i*i) + a0;

The resulting values to laster blaster's magic abs-condition are
c = r*(r2-3*i2);
d = a*(3*r2 + a2) + 3*r*(a2 - 2*i*b - b2) - 3*a*(i2 + 2*i*b + b2);

y-axis:
zi = abs(((zrsqr * 3.0) - zisqr) * zi) + JuliaI;
abs(((x*x * 3.0) - y*y) * y) + y0 - abs(((r*r * 3.0) - i*i) * i) - i0;
abs((3*r*r + 6*r*a + 3*a*a - i*i - 2*i*b - b*b) * (i+b)) + b0 - abs(3*r*r*i - i*i*i);
abs(3*r*r*i-i*i*i + 6*r*a*i + 3*a*a*i - 3*i*i*b - 3*i*b*b + 3*r*r*b + 6*r*a*b + 3*a*a*b - b*b*b) + b0 - abs(3*r*r*i-i*i*i);

And the resulting values to the magic abs-condition are:
c = i*(3*r2-i2);
d = 3*i*(2*r*a + a2 - b2) + 3*b*(r2 + 2*r*a + a2) - b*(3*i2 + b2);

Hey Kalles, thanks for the new types! Since today I'm exploring the Celtic type. Already a nice image (just a little after touch with the Gimp):

(http://fc02.deviantart.net/fs71/f/2014/280/6/2/celtic_curls_on_the_horizon_by_jeroensnake-d81ycqa.jpg)
Source: http://jeroensnake.deviantart.com/art/Celtic-Curls-on-the-horizon-487011826

Expect zoom videos the next few weeks!

:smileylineup:

Excellent, thanks!

Reposting here, because it needs to be repeated:

This is awesome news; thank you so much!!!

Here are the ABS fractal formulas I have cataloged so far. The 2nd order fractal formulas (set of 12 formulas) are complete:
(http://fc06.deviantart.net/fs70/f/2014/280/0/8/mandelbrot_abs_variations__complete_by_stardust4ever-d81yycs.png)
http://sta.sh/0ag7nkkbjri

3rd order has a long way to go (bigger superset but lots of duplicates and "junk") but here are some of the more interesting fractals I've discovered so far.
(http://fc07.deviantart.net/fs70/f/2014/280/b/0/cubic_mandelbrot_abs_variations__incomplete_by_stardust4ever-d81yyfd.png)
http://sta.sh/0263mwwvgaxn

Click the download link for full size view (1920x1920). Floating point code for each fractal has been embedded in the image.

I thought I ought to make the formulas public so that other people can play with them! Also I don't quite have the programming skills to incorporate them into functional plugins... :wink:

Hi! I have implemented perturbation for Perpendicular burning ship for GPGPU.  Initially I take the code from Kalles Fraktaler without changing (as it possible):
But I was surprised...
Why such constructions
and multiple copies the same numbers...  I understand that Kalles Fraktaler's author is very good programmer and maybe I don't know some important things. I have done some algebraic reconstructions and write this:
The pictures absolutely identical. New code 3x faster (1551 seconds against 4620).

Awesome if these advances can be built into the next KF build. Will GPU rendering work well on an Nvidia Gforce 430? 96 Cuda cores. My GPU is pretty old as I never found a need to upgrade it. Will it still outperform my AMD Bulldozer 8150 @4.2Ghz?

Yes, I use GeForce 730 (Fermi). It is the same GF108 (96 cores) that GeForce 430. Your results can be similar (IMHO)  88)

That is awesome that you are able to generate the abs fractals 3 times faster!
No, I haven't had time to try different conditions in that way :)

Maybe, it works well only in GPUs, because they are sensitive to different if-else conditions int the same warp. Also, GPUs are sensitive to number of variables, because register file is limited and latency of global memory is much more.

This is the code in my next version of Kalles Fraktaler (2.12.2, coming soon):


(and similarly for long double and floatexp)

One test is 7.5% faster than the original code in my fork, with no differences in the image.  A nice improvement, even though not as impressive as the difference on GPU.  Thanks therror!

Note that the "2.0 *" formulation is important, trying to use "cd + c" to save a multiplication leads to unacceptable rounding errors (different images, actually slower due to more references needed, just wrong I guess).