fractalforums.com - Non-Trigonometric Expansions for Cosine Formula

bugman

Iterator

Posts: 122

Non-Trigonometric Expansions for Cosine Formula

« on: November 23, 2009, 09:40:55 AM »

You can see see the non-trigonometric expansions for the "sine" formula here:
http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/msg8680/#msg8680

I still think the "sine" formula makes the most sense because it gives {x, y, z}^0 = {1, 0, 0}.

However, if I understand correctly, Daniel White and Garth Thornton have been experimenting with a "cosine" formula. The non-trigonometric expansions for the "cosine" formula are as follows:

CosineFormulas.gif (14.73 KB, 955x523 - viewed 626 times.)
« Last Edit: November 23, 2009, 08:50:29 PM by bugman »	Logged

xenodreambuie

Conqueror

Posts: 57

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #1 on: November 23, 2009, 11:27:03 AM »

Thanks, Paul. One correction for the cosine formula: phi=n.arccos(z/r). It doesn't affect the rest.

I like the sine formula for the most part because it has the 2D Julia or Mandelbrot set in the Z=0 plane, and this works especially well for power -2. I stumbled on the cosine option initially because it didn't require any parity tricks to work in inverse iteration with the trig formulas, and I like the power 2 Julias it creates. It's worth reiterating that there are two variants of the cosine; one with correct roots and one with some flipped roots, depending on whether one uses forward or inverse iteration, and whether parity is checked. The one I prefer aesthetically has some flipped roots. For higher powers the sine and cosine Julias look increasingly similar, except for the symmetry of the bulbs. For even powers the cosine has rows of bulbs aligned with longitudinal lines, while the sine formula has them alternating.


	Logged

Regards, Garth
http://xenodream.com

bugman

Iterator

Posts: 122

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #2 on: November 23, 2009, 08:52:19 PM »

Thank you for drawing that to my attention Garth. I have updated my previous post accordingly. The formula makes more sense now.


	Logged

twinbee

Fractal Phenom

Posts: 383

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #3 on: August 12, 2010, 04:51:49 PM »

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;
}


« Last Edit: August 12, 2010, 06:33:49 PM by twinbee »	Logged

Jesse

Download Section
Fractal Schemer

Posts: 995

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #4 on: August 12, 2010, 06:46:10 PM »

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 wink


« Last Edit: August 12, 2010, 07:51:11 PM by Jesse »	Logged

twinbee

Fractal Phenom

Posts: 383

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #5 on: August 13, 2010, 12:36:55 AM »

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?


« Last Edit: August 13, 2010, 12:45:09 AM by twinbee »	Logged

Jesse

Download Section
Fractal Schemer

Posts: 995

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #6 on: August 13, 2010, 07:07:19 PM »

Quote from: twinbee on August 13, 2010, 12:36:55 AM

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.

Quote

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

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
}


« Last Edit: August 14, 2010, 04:27:20 PM by Jesse »	Logged

Jesse

Download Section
Fractal Schemer

Posts: 995

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #7 on: August 14, 2010, 04:31:23 PM »

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
}


	Logged

twinbee

Fractal Phenom

Posts: 383

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #8 on: August 15, 2010, 10:29:57 PM »

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


	Logged

Jesse

Download Section
Fractal Schemer

Posts: 995

Re: Non-Trigonometric Expansions for Cosine Formula

« Reply #9 on: August 16, 2010, 12:01:17 AM »

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!


	Logged

Pages: [1] Go Down

« previous next »

	Author	Topic: Non-Trigonometric Expansions for Cosine Formula (Read 2654 times)
		Description:
0 Members and 1 Guest are viewing this topic.

Related Topics
	Subject	Started by	Replies	Views	Last post
	Secant method for cosine function Images Showcase (Rate My Fractal)	HPDZ	1	639	February 03, 2009, 09:51:16 PM by cKleinhuis
	Secant Cosine 2D Art	HPDZ	0	395	May 01, 2009, 01:32:06 AM by HPDZ
	New 3D IFS formula for UF5 IFS - Iterated Function Systems « 1 2 »	fractalrebel	18	13163	November 13, 2011, 01:39:30 AM by xenodreambuie
	The Formula Fractal Movies	subblue	0	490	March 30, 2010, 08:22:26 PM by subblue
	a formula for everything Non-Fractal related Chit-Chat « 1 2 »	teamfresh	16	1754	August 30, 2010, 07:25:50 PM by hermann

	Welcome, Guest. Please login or register., Guest. Please login or register.	May 23, 2013, 03:29:06 AM
		Login with username, password and session length