Logo by Trifox - Contribute your own Logo!
News: Visit the official fractalforums.com Youtube Channel
 
*
Welcome, Guest. Please login or register. October 22, 2014, 05:24:17 PM


Login with username, password and session length



Pages: [1]   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: Non-Trigonometric Expansions for Cosine Formula  (Read 3295 times)
0 Members and 1 Guest are viewing this topic.
bugman
Conqueror
*******
Posts: 122



WWW
« 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 748 times.)
« Last Edit: November 23, 2009, 08:50:29 PM by bugman » Logged
xenodreambuie
Safarist
******
Posts: 80



WWW
« 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
Conqueror
*******
Posts: 122



WWW
« 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 Fertilizer
*****
Posts: 383



WWW
« 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: 1018


« 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  smiley
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 Fertilizer
*****
Posts: 383



WWW
« 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: 1018


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

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

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: 1018


« 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 Fertilizer
*****
Posts: 383



WWW
« 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: 1018


« 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
  Print  
 
Jump to:  


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

Powered by MySQL Powered by PHP Powered by SMF 1.1.20 | SMF © 2013, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.401 seconds with 31 queries. (Pretty URLs adds 0.023s, 2q)