QUICK Octal and Hex Fractals

[Charleswehner]

It may seem obvious when pointed out, but ZZ⁸ is slower than Z((Z²)²)².

Similarly, ZZ¹⁶ is four times slower than Z(((Z²)²)²)².

I could never resist a bargain, so I modified my original Mandelbrot files to work this way. Lots of multiplications for the price of few.

Here is the octal set:

This was made with http://wehner.org/tools/fractals/oct/oct.asm

Paul N. Lee (Nahee Enterprises) has pointed out elsewhere in these forums
   http://www.fractalforums.com/index.php?topic=285.msg1317#msg1317
that there is always one segment missing as compared with the power of Z. So the eightfold power gives a seven-sided figure.

I moved into the image to 158,285 (pixels in from left and from bottom), and enlarged 64 times linear. This produced this:

This is perfectly standard. It is all part of the well-trodden path of classical Julia/Mandelbrot chaos mathematics - which I have classified as exponential exponentiation. The source is at http://wehner.org/tools/fractals/oct/oct2.asm

I home in on pixel 261,345 of this image to get a fixed vector for a Julia set. Here it is:

The source code is at http://wehner.org/tools/fractals/oct/joct2.asm

Now I zoom out by means of http://wehner.org/tools/fractals/oct/joct.asm:

Surprise, surprise - there are eight "petals" to the Julia set, whilst only seven in the Mandelbrot.

So by covering old ground, I was not wasting my time. I learned something I had not known.

For completeness, I made the hexadecimal Mandelbrot by means of http://wehner.org/tools/fractals/hex/hex.asm

There are fifteen segments.

I move in to 146,290 at 64 by 64 times enlargement, with http://wehner.org/tools/fractals/hex/hex2.asm

I select pixel 212,316 for http://wehner.org/tools/fractals/hex/jhex2.asm

I zoom out of the 64 times enlargement with http://wehner.org/tools/fractals/hex/jhex.asm

SIXTEEN segments.

Charles

[gandreas]

It may seem obvious when pointed out, but Z⁸ is slower than Z((Z²)²)².

It all depends on your math library - if you look, most "integer power" routines, rather than doing:

Code:

result = 1
for i = 1 to N
  result = result * Z

they actually go through the bits of N and do the far more optimal version:

Code:

if N is odd then result = Z else result = 1
N = N >> 1
while N > 0
  Z = Z * Z
  if N is odd then result = result * Z
  N = N >> 1

Floating point exponentiation pretty much requires using log, however...

Surprise, surprise - there are eight "petals" to the Julia set, whilst only seven in the Mandelbrot.

If you think about, it's just like a "regular" Mandelbrot, which only has one "petal" while the Julia Set has two for z².  Since you've been able to show how that first petal grows on the Mandelbrot set for functions that grow less, it might be interesting to do the same experiment plotting the Julia set (in general it seems that Julia sets can produce more interesting images for the same formula would be boring as a standard Mandelbrot style function).

Taking it one step further, you can view both as a single 4D objects (where the four coordinates are {Zr,Zi,Cr,Ci}) - the Mandelbrot set is then just a 2D cross section along the {Zr,Zi} axis, and the Julia set is a cross section along the {Cr,Ci} axis.  (The Budhabrot is then projection across the {Zr,Zi} axis which instead of being a cross section stochastically collapses the {Cr,Ci} axis).  Of course, you can get some interesting results by picking other axis pairs, or using various projections to rotate the 2D cross section across this shape for a more general case.

It would be fun to see which axis/projection can produce an "interesting" image using the "lowest" function.

[David Makin]

Floating point exponentiation pretty much requires using log, however...

Correct (for non-integer powers) but I'm pretty sure that in formulas where the power/degree in a calculation is a constant and an integer power of 2 then Frederik has written the UF compiler so that it uses the optimised method (using the FPU). In fact I'll bet he's probably done the compiler such that all constant integer powers up to say 15 or 31 are done using the optimum method (much further and its quicker to use the log method in any case).

[gandreas]

In fact I'll bet he's probably done the compiler such that all constant integer powers up to say 15 or 31 are done using the optimum method (much further and its quicker to use the log method in any case).

Not so sure about this - if, for example, we pick a simple "worse case" 255 (which, depending on your initial value and the size of the your floating point exponent, can easily overflow), this can be done with 7 floating point multiplications and 8 floating point addition:

Code:

total = z; // 1
z = z * z; 
total += z; // 3
z = z * z; 
total += z; // 7
z = z * z; 
total += z; // 15
z = z * z; 
total += z; // 31
z = z * z; 
total += z; // 63
z = z * z; 
total += z; // 127
z = z * z; 
total += z; // 255

While there are some funky ways to do log/exp with fixed point, unless you've got hardware log/exp instruction I don't believe you can beat this...

Of course, most of the time with fractals you end up doing complex math, so "z * z" is actually a complex multiply (since this algorithm works with anything where powers are defined as repeated multiplication), which is more expensive (four multiplies and two adds), but then again, complex log/exp are even more expensive as well...

[David Makin]

I was thinking of only FPU calculations - but you're probably correct it may be more optimum using multiply/add up to further than 31 even using the FPU (I haven't checked the relevant cycle times but I'm sure Frederik has) - obviously many special cases will always be faster using the multiply/add method - it depends on the complexity of the layout of the bits in the integer power.
When it comes to raising quaternions or hypercomplex numbers to integer powers, the multiply/add method becomes less optimum compared to using the log method.

[gandreas]

quaternions, in the general case, don't have logs or exp - only unit quaternions do.

So q ^ n = exp(log(q) * n) iff |q| == 1

Which means that the "optimized unrolled multiplication" is the only way to get integer powers of quaternions for all quaternions (and non-integer powers, well, given the non-commutative nature of quaternions, I'm not sure that they exist either).

[eNZedBlue]

To do fast exponentiation of complex numbers, I convert to polar co-ordinates and then multiply the polar angle by the exponentiation power, raise the vector length by the power, and then convert back to Cartesian co-ords. This relies on the geometric property of complex multiplication of two numbers on the complex plane being equivalent to adding the angles and multiplying the lengths of their polar representation.

Here is some simple C code for computing Z_n+1 = Z_n^power + C

C = dBi + dA
Z = dYi + dX
power = dPower

      double dU = dX * dX;
      double dV = dY * dY;
      double dLength = sqrt(dU + dV);
      double dAngle = atan2(dY, dX);

      dLength = pow(dLength, dPower);      
      dAngle = dAngle * dPower;

      dX = (cos(dAngle) * dLength) + dA;
      dY = (sin(dAngle) * dLength) + dB;

The nifty thing about this is that you can use arbitrary powers and it works in constant time. You can also do fractional powers, but if you do then you need to remember the angle from iteration to iteration, or it will get truncated in the transformation from polar to cartesian co-ordinates and back again (this also happens for integer powers, but is not an issue since the range always ends up being a multiple of pi radians). I'll post some videos of animating the exponent by fractional amounts once I get around to doing an over-night render. As you increase the exponent the set tends closer and closer to a circle, but detail is still evident if you zoom in.

Some renders:

Z_n+1 = Z_n¹⁰ + C (9 segments)

http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf10.jpg

Z_n+1 = Z_n⁵⁰⁰ + C (499 segments )

http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500.jpg

Z_n+1 = Z_n⁵⁰⁰ + C (detail)

http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500Detail.jpg

[lycium]

chris, unfortunately accurate division/sin/cos/pow/arctan are out of the question when using fixed point arithmetic - ironically this decision was made for "speed" reasons!

btw, you can save a square root in your code:

[r*e^(i*theta)]^n = r^n * e^(i*theta*n)

which boils down to

const double mod = pow(a*a + b*b,n * 0.5);
const double arg = atan2(b,a) * n;

const double new_a = mod * cos(arg);
const double new_b = mod * sin(arg);

[Nahee_Enterprises]

Chris "eNZedBlue" wrote:
>
>    Z_n+1 = Z_n⁵⁰⁰ + C (detail of the 499 segments )
>    http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500Detail.jpg


Interesting zoom into this area.

[Charleswehner]

It may seem obvious when pointed out, but Z⁸ is slower than Z((Z²)²)².

It all depends on your math library - if you look, most "integer power" routines, rather than doing:

Code:

result = 1
for i = 1 to N
  result = result * Z

&c.

That is why I do not use a math library. I simply code my own.

The internal prioritising of compilers and interpreters is not always published, so it takes time to familiarise oneself with other peoples' products. I prefer to move in, get the results and move on, with no fuss.

The nested squaring makes all power-of-two accessible. So I would not have gone for Z⁵⁰⁰ but rather Z⁵¹², for example - in the example from eNZedBlue:

http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500Detail.jpg

Twelve more powers may make a difference, but it is unlikely to be much.

Another concept, of course, is to make the exponentiation into a pastiche of the multiplication. For numbers up to 255

Register AH has the power, Register AL has nil.
Store Z
Shift right AX
Multiply Z by Z and store
If AH=0 then escape
Shift right AX
Multiply Z² by Z² and store
If AH=0 then escape
.........


After escaping, one will have the highest binary power in the registers, and all the half-powers and quarter-powers as well

Now one adds AL to itself (or shifts left). The top bit is lost. Thereafter:

Shift left AL
If carry, then multiply result by previous stored result
If AL is zero, then escape (DONE)
Shift left AL
.........


These are loops, of course. I opened them out to show the sequence of events. There would be seven repeats in the first loop, giving Z⁰ up to Z¹²⁸. Thereafter, there might be seven further multiplications - making a total of fourteen.

So with a maximum of fourteen multiplications, one can raise a complex number to the power of anything up to 255.

To do fast exponentiation of complex numbers, I convert to polar co-ordinates and then multiply the polar angle by the exponentiation power, raise the vector length by the power, and then convert back to Cartesian co-ords. This relies on the geometric property of complex multiplication of two numbers on the complex plane being equivalent to adding the angles and multiplying the lengths of their polar representation.

I read here on these forums that the Glynn set is Z to the power of one-and-a-half plus the vector.

If it is so, it is tricky. There are TWO square-roots, so one needs the correct one when multiplying by Z. Seen as a de Moivre polar manoevre, the angle should be added, giving one-and-a-half. If one makes a mistake, one subtracts and gets Z to the power a half.

Working with polar co-ordinates avoids such mistakes. Indeed, it would be quicker but for the VECTOR. Given ZZ^q, where q is some quantity, we just multiply angle by q and raise the range to the power of q. Given Given ZZ^q + Z_vector, we are in trouble. It is very tedious to convert back to Cartesian, using arcsine, arccos, arctan or whatever, and using the logarithm.

There is also - for the logarithm - the Problem of Multiple Solutions to consider.

I suspect that the Glynn set might be ZABS(Z)^3/2 + Z_vector.

Chris "eNZedBlue" wrote:

>    Z_n+1 = Z_n⁵⁰⁰ + C (detail of the 499 segments )
>    http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500Detail.jpg

Interesting zoom into this area.

Yes - even though we have Z to the power of 500 to the power of n, it is still in the family , which is exponential exponentiation. Z is exponentiated by m (which is 500), in turn exponentiated by n, the palette count.

Sunset on the Fractals seems to appear when one tries Mandelbrot for the first iteration (Z²), tribrot for the next (Z³) and tetrabrot, pentabrot and on - giving an initial term going Z², Z⁶, Z²⁴. That is the FACTORIAL exponentiation I explored elsewhere.

It came as quite a surprise to me, though, that this image (Z to the 500) is not more of the old, but introduces shapes not seen in the lower sets.

Charles

[lycium]

The nested squaring makes all power-of-two accessible. So I would not have gone for Z⁵⁰⁰ but rather Z⁵¹², for example - in the example from eNZedBlue:

http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500Detail.jpg

Twelve more powers may make a difference, but it is unlikely to be much.

in theory, but in practice (complex) multiplying with fixed point numbers 9 times will make the image really blocky; what's more it won't be anywhere near as fast the direct computation (to say nothing of the relative programming ease).

edit: btw,

There is also - for the logarithm - the Problem of Multiple Solutions to consider.

I suspect that the Glynn set might be Z<Quoted Image Removed>ABS(Z)^3/2 + Z_vector.

since all non-integer complex exponents are computed in terms of the logarithm, they will have infinitely many solutions (unique up to the choice of logarithm).

also: doing two square roots instead of one pow with 0.25 is much slower and inaccurate.

[Charleswehner]

I was sufficiently impressed by the image from Chris (eNZedBlue) to decide to do one of my own. This is what he showed us:

http://i19.photobucket.com/albums/b199/Stormfronter/MandelbrotPowerOf500Detail.jpg

Using nested squarings, we get 511 multiplications for the price of nine - a bargain - and don't even have to wait for the January sales. I modified my files to http://wehner.org/tools/fractals/512/512.asm and made this:

There are 511 teeth on this Mandelbrot cogwheel.

Then, with http://wehner.org/tools/fractals/512/512a.asm , I moved in closer and got this:

That image was from 144 pixels in and 320 up, and subsequently 64 fold enlarged.

Now I selected pixel 301,215 for a Julia vector and made http://wehner.org/tools/fractals/512/j512a.asm which gave this:

I've CRASHED THE GEARS! I never do that when out driving! It can't be me. It seems that the chaos mathematics is trying to fit 512 Julia teeth to a 511-tooth Mandelbrot hub.

The only tooth that has not been stripped, I conjecture, is that on the extreme left.

Charles

[David Makin]

quaternions, in the general case, don't have logs or exp - only unit quaternions do.

So q ^ n = exp(log(q) * n) iff |q| == 1

Which means that the "optimized unrolled multiplication" is the only way to get integer powers of quaternions for all quaternions (and non-integer powers, well, given the non-commutative nature of quaternions, I'm not sure that they exist either).

You are correct about the non-commutative problem but that simply means there's more than one possibility and in any case it applies equally to the multiply/add method as it does to the log/exp method :-)
For instance there are two ways of calculating sqrt(4) but that doesn't invalidate the math !

If I remember correctly I found that if the power is complex rather than fully quaternionic, it didn't matter anyway, certainly it doesn't matter if the power is real.

Here's an animation of the Newton for the roots of 1, from the cubic roots to the 30th roots:
http://www.deviantart.com/deviation/44441935/

Here's the code for q^n+c and the Newton q^n-r from my "Solid3D Quaternions" formula for Ultrafractal:

            elseif @fractal==5
              if (xi = sqrt(y*y + zz*zz + w*w))<1e-200
                x = x^real(@pwr) + @xconst
                y = @yconst
                zz = @zconst
                w = @wconst
              else
                z = log(x+flip(xi))
                zi=imag(z)/xi
                xi=real(z)
                yi = zi*y
                wi = zi*w
                zi = zi*zz
                x=real(@pwr)*xi-imag(@pwr)*yi-real(@pwr1)*zi-imag(@pwr1)*wi
                y=imag(@pwr)*xi+real(@pwr)*yi+imag(@pwr1)*zi-real(@pwr1)*wi
                zz=real(@pwr1)*xi-imag(@pwr1)*yi+real(@pwr)*zi+imag(@pwr)*wi
                w=imag(@pwr1)*xi+real(@pwr1)*yi-imag(@pwr)*zi+real(@pwr)*wi
                if (xi = sqrt(y*y + zz*zz + w*w))<1e-200
                  x = exp(x) + @xconst
                  y = @yconst
                  zz = @zconst
                  w = @wconst
                else
                  z = exp(x+flip(xi))
                  zi = imag(z)/xi
                  x = real(z) + @xconst
                  y = zi*y + @yconst
                  zz = zi*zz + @zconst
                  w = zi*w + @wconst
                endif
              endif
            elseif @fractal==6
              lx=x
              ly=y
              lz=zz
              lw=w
              if (xj = sqrt(y*y + zz*zz + w*w))<1e-200
                xi = x^real(pwr)
                yi = 0.0
                zi = 0.0
                wi = 0.0
              else
                z = log(x+flip(xj))
                zj = imag(z)/xj
                xj = real(z)
                yj = zj*y
                wj = zj*w
                zj = zj*zz
                xi=real(pwr)*xj-imag(pwr)*yj-real(@pwr1)*zj-imag(@pwr1)*wj
                yi=imag(pwr)*xj+real(pwr)*yj+imag(@pwr1)*zj-real(@pwr1)*wj
                zi=real(@pwr1)*xj-imag(@pwr1)*yj+real(pwr)*zj+imag(pwr)*wj
                wi=imag(@pwr1)*xj+real(@pwr1)*yj-imag(pwr)*zj+real(pwr)*wj
                if (xj = sqrt(yi*yi + zi*zi + wi*wi))<1e-200
                  xi = exp(xi)
                  yi = 0.0
                  zi = 0.0
                  wi = 0.0
                else
                  z = exp(xi+flip(xj))
                  zj = imag(z)/xj
                  xi = real(z)
                  yi = zj*yi
                  zi = zj*zi
                  wi = zj*wi
                endif
              endif
              xk=x*xi-y*yi-zz*zi-w*wi
              yk=y*xi+x*yi+w*zi-zz*wi
              zk=zz*xi-w*yi+x*zi+y*wi
              wk=w*xi+zz*yi-y*zi+x*wi
              if @mulfirst
                xj=real(@pwr)*xi-imag(@pwr)*yi-real(@pwr1)*zi-imag(@pwr1)*wi
                yj=imag(@pwr)*xi+real(@pwr)*yi+imag(@pwr1)*zi-real(@pwr1)*wi
                zj=real(@pwr1)*xi-imag(@pwr1)*yi+real(@pwr)*zi+imag(@pwr)*wi
                wj=imag(@pwr1)*xi+real(@pwr1)*yi-imag(@pwr)*zi+real(@pwr)*wi
                xi=real(pwr)*xk-imag(pwr)*yk-real(@pwr1)*zk-imag(@pwr1)*wk \
                   +@xconst
                yi=imag(pwr)*xk+real(pwr)*yk+imag(@pwr1)*zk-real(@pwr1)*wk \
                   +@yconst
                zi=real(@pwr1)*xk-imag(@pwr1)*yk+real(pwr)*zk+imag(pwr)*wk \
                   +@zconst
                wi=imag(@pwr1)*xk+real(@pwr1)*yk-imag(pwr)*zk+real(pwr)*wk \
                   +@wconst
              else
                xj=xi*real(@pwr)-yi*imag(@pwr)-zi*real(@pwr1)-wi*imag(@pwr1)
                yj=yi*real(@pwr)+xi*imag(@pwr)+wi*real(@pwr1)-zi*imag(@pwr1)
                zj=zi*real(@pwr)-wi*imag(@pwr)+xi*real(@pwr1)+yi*imag(@pwr1)
                wj=wi*real(@pwr)+zi*imag(@pwr)-yi*real(@pwr1)+xi*imag(@pwr1)
                xi=xk*real(pwr)-yk*imag(pwr)-zk*real(@pwr1)-wk*imag(@pwr1) \
                   +@xconst
                yi=yk*real(pwr)+xk*imag(pwr)+wk*real(@pwr1)-zk*imag(@pwr1) \
                   +@yconst
                zi=zk*real(pwr)-wk*imag(pwr)+xk*real(@pwr1)+yk*imag(@pwr1) \
                   +@zconst
                wi=wk*real(pwr)+zk*imag(pwr)-yk*real(@pwr1)+xk*imag(@pwr1) \
                   +@wconst
              endif
              sq=1.0/(xj*xj+yj*yj+zj*zj+wj*wj)
              x=(xi*xj+yi*yj+zi*zj+wi*wj)*sq
              y=(yi*xj-xi*yj-wi*zj+zi*wj)*sq
              zz=(zi*xj+wi*yj-xi*zj-yi*wj)*sq
              w=(wi*xj-zi*yj+yi*zj-xi*wj)*sq

[dentaku2]

So, if we talk about speed, how long did your images take? (My generator: (root Mandelbrot, 800x600, 100 iterations, exponent = 500): 28 seconds.)

[lycium]

So, if we talk about speed, how long did your images take? (My generator: (root Mandelbrot, 800x600, 100 iterations, exponent = 500): 28 seconds.)

two points to make here:

• experience tells that you should really not mention speed in this thread (even if it's about assembly language and integers)
• to give a benchmark without so much as a machine spec is wholly useless