Triplex algebra

[David Makin]

Multiplication on cartesian coordinates is not associative or distributive *because first and foremost it makes no sense*.  Multiplication on spherical coordinates (i.e. multiply radius, sum angles) or triplex-expressed-as-quaternion is commutative and associative, not sure about distributive.

I'm afraid you've lost me there.

For the triplex maths as applied in the Mandelbulb, the following are *not* true for all values:

a*(b*c) = a*b*c = (a*b)*c

a*(b+c) = a*b + a*c

a*(b/a) = b

exp(ln(a)) = a

If for your quaternion form any of these *do* hold true then it's not the same mathematical form as that used for the original Mandelbulb - I'm not saying it can't produce similar fractals though.

Note that I think however that the above equalities obviously do hold true if in a, b and c the "j" component is zero i.e. we're reduced to complex numbers.

I should add that even:

(a^2)*a = a^3

does not hold true for all values.

[David Makin]

Multiplication on cartesian coordinates is not associative or distributive *because first and foremost it makes no sense*.  Multiplication on spherical coordinates (i.e. multiply radius, sum angles) or triplex-expressed-as-quaternion is commutative and associative, not sure about distributive.

I'm afraid you've lost me there.

For the triplex maths as applied in the Mandelbulb, the following are *not* true for all values:

a*(b*c) = a*b*c = (a*b)*c

a*(b+c) = a*b + a*c

a*(b/a) = b

exp(ln(a)) = a

If for your quaternion form any of these *do* hold true then it's not the same mathematical form as that used for the original Mandelbulb - I'm not saying it can't produce similar fractals though.

Note that I think however that the above equalities obviously do hold true if in a, b and c the "j" component is zero i.e. we're reduced to complex numbers.

I should add that even:

(a^2)*a = a^3

does not hold true for all values.

Just been thinking about it and realised that the Paolo is correct with respect to treating the triplex in polar terms when dealing solely with multiplication/division since if you convert all the terms to polar first and then perform the calculations then associativity works and a*(b/a)=a *but* as soon as you need to include addition/subtraction then things break down because however you do the calculations (in 3 dimensions only) the distributive law does not apply.

I'm just wondering if the derivations done in this thread for the higher functions using Mathematica etc. were done in such a way that Mathematica took into account that the distributive law does not apply ?
Or were they all done using Paolo's quaternion form ?

An obvious problem is something as basic as the derivation of the derivative of y = x^2 where:

y+dy = (x+dx)^2
y + dy = x^2 + 2*x*dx + dx^2
x^2 + dy = x^2 +2*x*dx + dx^2
dy = 2*x*dx + dx^2
dy/dx = 2*x + dx
So as dx->0:
Dy/Dx = 2*x

The problem being that for the triplex you cannot expand (x+dx)^2 as x^2 + 2*x*dx + dx^2 for all x/dx.

This may be part of the reason why the analytical DE seems to break down in certain areas.

[Timeroot]

Okay, so it *is* associative, since you're keeping it in spherical coordinates and only using addition/multiplication. Also, the division/multiplication inverse *does* work, agreed? I can understand, though, why the distributive property doesn't work. I'm curious as to how exactly exp(a) is defined - you could use taylor series, or the exponentiation formula I suggested with A=2.71..., or exp(x)=e^x * (cos x,sin y,0) * (cos z, 0, sin z) [or WHATEVER that was on the first page].... it seems there is an agreed upon definition for Mandelbulb purposes, can someone please tell me what it is?

With regard to derivatives, I think it would be interesting how exactly (x+dx)^2 can be expanded [that is, if it can be].

EDIT:
 Forgive me if I'm being a pain with having all this explained, but why doesn't (z*z)*z=z^3? I'd think that z^3=Result, where Result.Radius=z.Radius^3, Result.Theta=z.Theta*3, and Result.Phi=z.Phi*3. z*z=ResultA, where ResultA.Radius=z.Radius*z.Radius, ResultA.Theta=z.Theta + z.Theta, and and ResultA.Phi=z.Phi + z.Phi. Then (z*z)*z=ResultA*z=Result B, where ResultB.Radius=z.Radius * ResultA.Radius = (substituting here) z.Radius * z.Radius * z.Radius = z.Radius^3. ResultB.Theta=z.Theta + ResultA.Theta = z.Theta + z.Theta + z.Theta = z.Theta*3. Same goes for Result.Phi

Clearly, z^3=Result should equal (z*z)*z=ResultB. Can somebody explain to me why it is not?

[David Makin]

Okay, so it *is* associative, since you're keeping it in spherical coordinates and only using addition/multiplication. Also, the division/multiplication inverse *does* work, agreed? I can understand, though, why the distributive property doesn't work. I'm curious as to how exactly exp(a) is defined - you could use taylor series, or the exponentiation formula I suggested with A=2.71..., or exp(x)=e^x * (cos x,sin y,0) * (cos z, 0, sin z) [or WHATEVER that was on the first page].... it seems there is an agreed upon definition for Mandelbulb purposes, can someone please tell me what it is?

With regard to derivatives, I think it would be interesting how exactly (x+dx)^2 can be expanded [that is, if it can be].

EDIT:
 Forgive me if I'm being a pain with having all this explained, but why doesn't (z*z)*z=z^3? I'd think that z^3=Result, where Result.Radius=z.Radius^3, Result.Theta=z.Theta*3, and Result.Phi=z.Phi*3. z*z=ResultA, where ResultA.Radius=z.Radius*z.Radius, ResultA.Theta=z.Theta + z.Theta, and and ResultA.Phi=z.Phi + z.Phi. Then (z*z)*z=ResultA*z=Result B, where ResultB.Radius=z.Radius * ResultA.Radius = (substituting here) z.Radius * z.Radius * z.Radius = z.Radius^3. ResultB.Theta=z.Theta + ResultA.Theta = z.Theta + z.Theta + z.Theta = z.Theta*3. Same goes for Result.Phi

Clearly, z^3=Result should equal (z*z)*z=ResultB. Can somebody explain to me why it is not?

I guess it's my fault for not following Paolo's comments well enough - I was describing the issue if you calculate z^2 getting the result in cartesian form then back to polar to multiply by the original z giving z^2*z in cartesian form at the end or using the fully cartesian calculation to get z^2 and then say polar and back for z^2*z.
These do not result in the same value as doing the calculation of z^3 in polar form and only converting to cartesian at the end (at least not for all z).

[Paolo Bonzini]

Okay, so it *is* associative, since you're keeping it in spherical coordinates and only using addition/multiplication.

Only using multiplication/division, actually.  Addition requires cartesian coordinates, so it is associative but not distributive.

Also, the division/multiplication inverse *does* work, agreed? I can understand, though, why the distributive property doesn't work. I'm curious as to how exactly exp(a) is defined - you could use taylor series, or the exponentiation formula I suggested with A=2.71..., or exp(x)=e^x * (cos x,sin y,0) * (cos z, 0, sin z) [or WHATEVER that was on the first page].... it seems there is an agreed upon definition for Mandelbulb purposes, can someone please tell me what it is?

I have no idea.  The whatever was on the first page probably does not agree with the Taylor series, which is a problem.  Maybe quaternions could help, I don't have much time now.

Forgive me if I'm being a pain with having all this explained, but why doesn't (z*z)*z=z^3?

When you do z*z*z*z*z in cartesian form, every multiplication includes an implicit conversion to and from spherical coordinates.  It's a bit hard to visualize, but the conversion loses information because of the different range of the elevation vs. the y-axis rotation term.  Instead, you have to use specially crafted formulas for each exponent(*) that include exactly one conversion to spherical coordinates and one from, which is what is called z^n.

(*)I had a generic exponentiation formula expressed using Chebyshev polynomials, but I haven't looked at it for a while and I'm not sure I had no calculation errors since I did it on paper.

[twinbee]

Thanks to Paul for these! I had another crack at getting this to work, and finally managed to get it work yesterday. Ages ago, I tried replacing all y variables inside the If section with zero (apart from the main If condition obviously). But it turns out that replacing them with the same VeryLow number will make it work.

Here below is the code I use. It's tricky to tell because of overhead, but I find the Mandelbulb function with the below non-trig code generally around twice as fast as the trig version despite how long-winded it seems. Replacing the repetitive multiplications with precalculated variables doesn't seem to speed it up at all. I get another 10%+ speed boost, if in the: "if( fabs(y) < 0.000001 )" section, I use instead: a.x=0; a.y=0; a.z=0;  though this seems to produce a line across the center if the rotation of the scene is zero. The fact that I can replace x,y,z with such trivial numbers and it produce such nearly perfect results makes me believe there's a better way than this.

Code:

point tpow8(double x, double y, double z) {	
	point p;

	if( fabs(y) < 0.000001 ) {
		double r=sqrt(x*x+0.000001);
		double d=8*z*(z*z*z*z*z*z - 7*z*z*z*z*r*r + 7*z*z*r*r*r*r - r*r*r*r*r*r) / (r*r*r*r*r*r*r);
		p.x = d*(x*x*x*x*x*x*x*x);
		p.y = 0.000001;
		p.z = z*z*z*z*z*z*z*z-28*z*z*z*z*z*z*r*r+70*z*z*z*z*r*r*r*r-28*z*z*r*r*r*r*r*r+r*r*r*r*r*r*r*r;
	} else {
		double r=sqrt(x*x+y*y);
		double d=8*z*(z*z*z*z*z*z - 7*z*z*z*z*r*r + 7*z*z*r*r*r*r - r*r*r*r*r*r) / (r*r*r*r*r*r*r);
		p.x = d*(x*x*x*x*x*x*x*x - 28*x*x*x*x*x*x*y*y + 70*x*x*x*x*y*y*y*y - 28*x*x*y*y*y*y*y*y + y*y*y*y*y*y*y*y);
		p.y = 8*d*x*y*(x*x*x*x*x*x - 7*x*x*x*x*y*y + 7*x*x*y*y*y*y - y*y*y*y*y*y);
		p.z = z*z*z*z*z*z*z*z - 28*z*z*z*z*z*z*r*r + 70*z*z*z*z*r*r*r*r - 28*z*z*r*r*r*r*r*r + r*r*r*r*r*r*r*r;
	}
	return p;
}

[Jesse]

hi twinbee,

i also made just the quadratic temporary values, multiplications are done very fast (about 4 clks) on the cpu.
Going through memory is often slower, and i did for the cosine and sine versions also pure x87 asm, to do it with sse2 was a bit to much work for me  
If you need fast code and can get use of asm, then i can send you also the cosine code, on this computer i have only the sine version.
The pow8 anniversary is upcoming..

Code:

asm                       //Sine pow8 bulb
  push  esi
  push  edi
  mov   esi, [ebp + 8]    //PIteration3D
  fld   qword [eax]       //x
  mov   edi, [esi + 48]   //PVars
  fmul  st(0), st(0)      //xx
  fld   qword [edx]       //y
  add   edi, 88
  fmul  st(0), st(0)      //yy,xx
  fld   qword [ecx]       //z,yy,xx
  fmul  st(0), st(0)      //zz,yy,xx
  fld   st(2)             //xx,zz,yy,xx
  fadd  st(0), st(2)      //xx+yy=r,zz,yy,xx
  fld   st(0)             //r,r,zz,yy,xx
  fmul  st(0), st(1)      //rr,r,zz,yy,xx
  fld   st(2)
  fmul  st(0), st(0)      //zzzz(S3*S3),rr,r,zz,yy,xx
  fld   st(2)             //r,zzzz(S3*S3),rr,r,zz,yy,xx    z calculation
  fmul  st(0), st(4)      //r*zz
  fmul  qword [edi + 56]  //6*r*zz,zzzz(S3*S3),rr,r,zz,yy,xx   
  fsubr st(0), st(1)      //zzzz-6rzz,zzzz,rr,r,zz,yy,xx
  fadd  st(0), st(2)      //zzzz-6rzz+rr,zzzz,rr,r,zz,yy,xx
  fld   st(4)             //zz,zzzz-6rzz+rr,zzzz,rr,r,zz,yy,xx
  fsub  st(0), st(4)      //zz-r,zzzz-6rzz+rr,zzzz,rr,r,zz,yy,xx
  fmulp                   //(zz-r)*(zzzz-6rzz+rr),zzzz,rr,r,zz,yy,xx
  fld   st(3)             //r,(zz-r)*(zzzz-6rzz+rr),zzzz,rr,r,zz,yy,xx
  fsqrt
  fmulp                   //sqrt(r)*(zz-r)*(zzzz-6rzz+rr),zzzz,rr,r,zz,yy,xx
  fmul  qword [ecx]       //*z
  fmul  qword [edi + 72]  //*8
  fmul  qword [edi - 104] //*dZmul
  fchs
  fadd  qword [esi + 40]  //+J3
  fstp  qword [ecx]       //zzzz,rr,r,zz,yy,xx
  fld   st(0)             //zzzz,zzzz,rr,r,zz,yy,xx        a calculation
  fadd  st(0), st(2)      //zzzz+rr,zzzz,rr,r,zz,yy,xx
  fmulp st(3), st(0)      //zzzz,rr,r*(zzzz+rr),zz,yy,xx
  fld   st(1)             //rr,zzzz,rr,r*(zzzz+rr),zz,yy,xx
  fmul  qword [edi + 120]  //rr*70,zzzz,rr,r*(zzzz+rr),zz,yy,xx
  fadd  st(0), st(1)
  fmulp                   //(rr*70+zzzz)*zzzz,rr,r*(zzzz+rr),zz,yy,xx
  fxch  st(2)             //r*(zzzz+rr),rr,(rr*70+zzzz)*zzzz,zz,yy,xx
  fmulp st(3), st(0)      //rr,(rr*70+zzzz)*zzzz,zz*r*(zzzz+rr),yy,xx
  fxch  st(2)             //zz*r*(zzzz+rr),(rr*70+zzzz)*zzzz,rr,yy,xx
  fmul  qword [edi + 104]  //28*zz*r*(zzzz+rr),(rr*70+zzzz)*zzzz,rr,yy,xx
  fsubp                   //(rr*70+zzzz)*zzzz-28*zz*r*(zzzz+rr),rr,yy,xx
  fxch  st(1)
  fmul  st(0), st(0)      //rrrr,(rr*70+zzzz)*zzzz-28*zz*r*(zzzz+rr),yy,xx
  fdivp                   //(zzzz*(rr*70+zzzz-28*zz*r*(zzzz+rr))/rrrr,yy,xx
  fadd  qword [edi - 56]  //a,yy,xx    +1
  fld   st(1)             //yy,a,yy,xx                   y calculation
  fmul  qword [edi + 64]  //7*yy,a,yy,xx      
  fld   st(3)             //xx,7*yy,a,yy,xx
  fmul  qword [edi + 64]  //7*xx,7*yy,a,yy,xx
  fsub  st(0), st(3)      //7*xx-yy,7*yy,a,yy,xx
  fld   st(4)             //xx,7*xx-yy,7*yy,a,yy,xx
  fsubr st(2), st(0)      //xx,7*xx-yy,xx-7*yy,a,yy,xx
  fmul  st(0), st(0)      //xxxx,7*xx-yy,xx-7*yy,a,yy,xx
  fmul  st(2), st(0)      //xxxx,7xx-yy,xxxx(xx-7yy),a,yy,xx
  fld   st(4)             //yy,xxxx,7xx-yy,xxxx(xx-7yy),a,yy,xx
  fmul  st(0), st(0)      //yyyy,xxxx,7xx-yy,xxxx(xx-7yy),a,yy,xx
  fmul  st(2), st(0)      //yyyy,xxxx,yyyy(7xx-yy),xxxx(xx-7yy),a,yy,xx
  fxch  st(2)             //yyyy(7xx-yy),xxxx,yyyy,xxxx(xx-7yy),a,yy,xx
  faddp st(3), st(0)      //xxxx,yyyy,yyyy(7xx-yy)+xxxx(xx-7yy),a,yy,xx
  fxch  st(2)             //yyyy(7xx-yy)+xxxx(xx-7yy),yyyy,xxxx,a,yy,xx
  fmul  qword [edi + 72]  //*8
  fmul  qword [eax]       //*x
  fmul  qword [edx]       //*y
  fmul  st(0), st(3)      //*a
  fadd  qword [esi + 32]  //+J2
  fstp  qword [edx]       //yyyy,xxxx,a,yy,xx
  fld   st(1)             //xxxx,yyyy,xxxx,a,yy,xx
  fmul  qword [edi + 120]  //70xxxx,yyyy,xxxx,a,yy,xx
  fadd  st(0), st(1)      //70xxxx+yyyy,yyyy,xxxx,a,yy,xx
  fmul  st(0), st(1)      //yyyy(70xxxx+yyyy),yyyy,xxxx,a,yy,xx
  fxch  st(1)             //yyyy,yyyy(70xxxx+yyyy),xxxx,a,yy,xx
  fadd  st(0), st(2)      //yyyy+xxxx,yyyy(70xxxx+yyyy),xxxx,a,yy,xx
  fmulp st(4), st(0)      //yyyy(70xxxx+yyyy),xxxx,a,yy(yyyy+xxxx),xx
  fxch  st(4)             //xx,xxxx,a,yy(yyyy+xxxx),yyyy(70xxxx+yyyy)
  fmulp st(3), st(0)      //xxxx,a,xxyy(yyyy+xxxx),yyyy(70xxxx+yyyy)
  fmul  st(0), st(0)      //xxxx*xxxx,a,xxyy(yyyy+xxxx),yyyy(70xxxx+yyyy)
  faddp st(3), st(0)      //a,xxyy(yyyy+xxxx),xxxx*xxxx+yyyy(70xxxx+yyyy)
  fxch  st(1)             //xxyy(yyyy+xxxx),a,xxxx*xxxx+yyyy(70xxxx+yyyy)
  fmul  qword [edi + 104]
  fsubp st(2), st(0)      //a,xxxx*xxxx+yyyy(70xxxx+yyyy)-28xxyy(yyyy+xxxx)
  fmulp
  fadd  qword [esi + 24]
  fstp  qword [eax]
  pop   edi
  pop   esi
end

forgot to mention that i dont do a div0 test at all, i masked the exceptions and slightly rotated the bulb on startup, so this is not really a problem for me  

[twinbee]

Hey that looks promising. Any idea how much faster the SS2 version would run?

I haven't started learning assembler yet (I know to inline with __asm {...} ). How can I communicate the C++ x,y,z variables to the asm code? And does the asm code produce nx,ny,nz variables by the end or something?

[Jesse]

Hey that looks promising. Any idea how much faster the SS2 version would run?

A wild guess would be 50%, that depends a lot of how the code can be optimized for sse2.
Cant predict it because i just tried a minute and then had, or wanted to have, other things to do.  

Wish i had a smart compiler that could do it for me.

I haven't started learning assembler yet (I know to inline with __asm {...} ). How can I communicate the C++ x,y,z variables to the asm code? And does the asm code produce nx,ny,nz variables by the end or something?

In my case the code is a whole procedure with the beginning:

 //P8 sine bulb       x:eax, y:edx, z:ecx, w:esp->ebp+12, PIt:ebp+8
procedure HybridIntP8(var x, y, z, w: Double; PIteration3D: TPIteration3D);
asm
..
end

You can also inline code and use variables and constants directly, for example with

  fadd  qword Cx

where Cx is a double const or var, instead of

  fadd  qword [esi + 24]

where esi is as example a pointer to an array of double values, and +24 selects the fourth value.

I had to do it the last way because of the independence of the code, the PIteration3D points to fixed structure.


If you also want to do a whole procedure, then you need to know some points:

The parameters are submitted by specific registers that are depending on the calling conventions of the operating system and/or compiler language, just know the one i use... but that can be estimated for example by looking at some assembled code.

For the procedure above, the inputs are 4 pointer to the double values plus a pointer to a structure with constants like Cx,Cy,Cz and also another pointer to formula dependend variables and constants.

Hint: it may be better to use only one pointer to a 4d vector, so you need only 1 register instead of 4.
Once tested it in an older program version that was not asm optimized and did not get a speed benefit.
Maybe with handwritten code it would be better, but i cannot change it all now.


OK, for first tests the inline variant would be easier i guess.  

Maybe this a good starting point, though c is not my main language...
try and pray

Code:

double Zx, Zy, Zz, Cx, Cy, Cz;
double c6 = 6.0, c7 = 7.0, c8 = 8.0, c28 = 28.0, c70 = 70.0, Cvsv = 1e-40, dZmul = -1.0;

__asm                       // Sine pow8 bulb
{ fld   qword Zx   
  fmul  st(0), st(0)     
  fld   qword Zy
  fmul  st(0), st(0)     
  fld   qword Zz
  fmul  st(0), st(0)     
  fld   st(2)             
  fadd  st(0), st(2)      
  fld   st(0)             
  fmul  st(0), st(1)      
  fld   st(2)
  fmul  st(0), st(0)     
  fld   st(2)             
  fmul  st(0), st(4)     
  fmul  qword C6          // 6.0 double constant 
  fsubr st(0), st(1)      
  fadd  st(0), st(2)      
  fld   st(4)             
  fsub  st(0), st(4)      
  fmulp                   
  fld   st(3)             
  fsqrt
  fmulp                   
  fmul  qword Zz 
  fmul  qword c8           // 8.0 double constant
  fmul  qword dZmul      // double var -1.0 or 1.0 to change between +Z and -Z variant
  fchs
  fadd  qword Cz          // double Cz constant
  fstp  qword Zz       
  fld   st(0)             
  fadd  st(0), st(2)     
  fmulp st(3), st(0)     
  fld   st(1)             
  fmul  qword c70         // 70.0 double constant
  fadd  st(0), st(1)
  fmulp                  
  fxch  st(2)            
  fmulp st(3), st(0)      
  fxch  st(2)            
  fmul  qword c28         // 28.0 double constant
  fsubp                   
  fxch  st(1)
  fmul  st(0), st(0)     
  fadd  qword Cvsv       // a very small value to avoid div0
  fdivp   
  fld1              
  faddp  
  fld   st(1)             
  fmul  qword c7           // 7.0 double constant     
  fld   st(3)             
  fmul  qword c7           // 7.0 double constant
  fsub  st(0), st(3)      
  fld   st(4)             
  fsubr st(2), st(0)     
  fmul  st(0), st(0)      
  fmul  st(2), st(0)      
  fld   st(4)            
  fmul  st(0), st(0)      
  fmul  st(2), st(0)      
  fxch  st(2)            
  faddp st(3), st(0)      
  fxch  st(2)            
  fmul  qword c8           // 8.0 double constant
  fmul  qword Zx         
  fmul  qword Zy      
  fmul  st(0), st(3)  
  fadd  qword Cy
  fstp  qword Zy         
  fld   st(1)             
  fmul  qword c70         // 70.0 double constant
  fadd  st(0), st(1)      
  fmul  st(0), st(1)      
  fxch  st(1)             
  fadd  st(0), st(2)     
  fmulp st(4), st(0)      
  fxch  st(4)            
  fmulp st(3), st(0)      
  fmul  st(0), st(0)     
  faddp st(3), st(0)     
  fxch  st(1)             
  fmul  qword c28         // 28.0 double constant
  fsubp st(2), st(0)      
  fmulp
  fadd  qword Cx
  fstp  qword Zx
}

[Jesse]

1e-20 was way to big, 1e-40 is better but it does not really prevent from overflows, i fear.
I changed this in the previous post.

Here is the cosine code, you may not have to write qword in front of the double values because the compiler does know the length.
But my compiler is not very smart, so i better tell him what to do.
This code is on my PC 43% faster then the pascal version...

Code:

__asm{                    // cosine pow8
  fld   qword Zx
  fmul  st(0), st(0)   
  fld   qword Zy
  fmul  st(0), st(0)     
  fld   qword Zz 
  fmul  st(0), st(0) 
  fld   st(2)        
  fadd  st(0), st(2) 
  fld   st(0)         
  fmul  st(0), st(0)  
  fld   st(2)                                 
  fmul  st(0), st(0)   
  fld   st(1)             //  z calculation
  fadd  st(0), st(1)     
  fmul  qword c28
  fmul  st(0), st(3)
  fmul  st(0), st(4)     
  fld   st(1)            
  fmul  st(0), st(0)
  fsubrp                  
  fld   st(1)            
  fmul  st(0), st(3)
  fmul  qword c70
  faddp                  
  fld   st(2)
  fmul  st(0), st(0)
  faddp                  
  fmul  qword dZmul              
  fadd  qword Cz
  fld   st(3)             //  a calculation
  fsub  st(0), st(5)     
  fmul  qword c7
  fmul  st(0), st(5)     
  fsub  st(0), st(3)
  fmul  st(0), st(4)    
  fxch  st(2)            
  fmulp st(5), st(0)     
  fxch  st(4)            
  faddp                   
  fmul  qword Zx
  fmul  qword c8
  fld   st(2)
  fsqrt                  
  fmulp st(2), st(0)    
  fxch  st(2)            
  fmulp                   
  fadd  qword Cvsv
  fdivp                  
  fxch  st(1)
  fstp  qword Zz       
  fld   st(1)             //  y calculation
  fmul  qword c7  
  fld   st(3)        
  fmul  qword c7
  fsub  st(0), st(3)   
  fld   st(4)          
  fsubr st(2), st(0)  
  fmul  st(0), st(0)  
  fmul  st(2), st(0) 
  fld   st(4)        
  fmul  st(0), st(0)     
  fmul  st(2), st(0)     
  fxch  st(2)            
  faddp st(3), st(0)     
  fxch  st(2)            
  fmul  qword c8
  fmul  qword Zx
  fmul  qword Zy
  fmul  st(0), st(3)    
  fadd  qword Cy
  fstp  qword Zy      
  fld   st(1)             //  x calculation
  fmul  qword c70  
  fadd  st(0), st(1)     
  fmul  st(0), st(1)     
  fxch  st(1)            
  fadd  st(0), st(2)      
  fmulp st(4), st(0)     
  fxch  st(4)            
  fmulp st(3), st(0)    
  fmul  st(0), st(0)     
  faddp st(3), st(0)     
  fxch  st(1)            
  fmul  qword c28
  fsubp st(2), st(0)     
  fmulp
  fadd  qword Cx
  fstp  qword Zx
}

[twinbee]

I wonder how much faster your asm code runs compared to the C version. Anyway, like you hinted at, I removed any "qword" commands as this gave errors like "error C2400: inline assembler syntax error in 'first operand'; found 'Zx'".

However, now I'm getting errors like "warning C4410: illegal size for operand" for commands like "fmul   C6".

I also got multiple errors for "fmulp" commands like: "error C2414: illegal number of operands".

I'm guessing this is because there are slightly different flavours of assembly language, and Visual C++ obviously uses its own.

No worries though if it's awkward to convert...

[Jesse]

Visual c++ is afaik using masm, maybe someone can correct the code for it!

I did not found very much helpful things on the net, so it could be something like
  fmul qword ptr [c6]

or c6 is not a valid name to choose, or dunno.

  fmulp
or other x87 instructions with omitted operands means
  fmulp st(0), st(1)

An operation on the first and second register.

But here must be several people who knows it!

[soler]

Bugman has defined a triplex polar form using a particular matrix product. I have extended this idea to give 48 different polar forms. If you are interested then please have a look at: http://soler7.com/Fractals/Matrices%20to%20Triplex.pdf

[M Benesi]

Something a hell of a lot easier is.. the complex triplex (from another post of mine elsewhere in the forums).  It really makes the mathematical relationship to the original 2d set apparent:

You don't need anything too complex** to do triplex "algebra".  You simply require:

A) a square root function to calculate a magnitude:
  r1= sqrt(y^2+z^2)

B) a complex power function:
  complex_1= (x + i r1)^n
  complex_2= (y + i z)^n

C) a real power function:
  r3=r1^-n       (you are applying the magnitude of y and z two times (once in each complex number), so need to divide it out once)

D) the ability to directly access the real and imaginary components of the complex numbers:
  new x = real part of complex 1   + x pixel value  OR  x Julia seed   (for Julias use the seed, Mandys use the pixel value)
  new y = imaginary part of complex_1 * real part of complex 2 * r3  + y pixel value  OR  y Julia seed
  new z = imaginary part of complex_1 * imaginary part of complex 2 * r3  + z pixel value  OR  z Julia seed

  It's easily extended to higher dimensions... and is faster than the trig version in certain compilers (although I haven't tried them all).

  ** pun was and is still intended.... 2 complex..  although I didn't mention that it was intentional in the original post, as I felt it was a bit heavy handed to point out the pun.  This, however, has changed. 

[soler]

Here is an image from one of the 48 variations:

Many more are at
http://soler7.com/Fractals/3D0.html
All were made using Visins of Chaos.