Most likely, you know the following ratios from school:

$$

$$\ sin (\ alpha + \ beta) = \ sin \ alpha \ times \ cos \ beta + \ cos \ alpha \ times \ sin \ beta \\ \ cos (\ alpha + \ beta) = \ cos \ alpha \ times \ cos \ beta - \ sin \ alpha \ times \ sin \ beta$$

When you first became acquainted with this formula in childhood, most likely, your first feeling was pain due to the fact that this formula must be remembered. This is very bad, because in fact you do not need to remember this formula - it is displayed when you rotate the triangle on paper. In fact, I do the same when I write down this formula. This interpretation will be apparent by the middle of this article. But now, to leave all the fun for later and to postpone the moment when you say “Eureka!”, Let's think about why we should even think about this formula.

Introduction

The trigonometric functions sin()



and cos()



probably the most popular in computer graphics, since they are the basis for describing any round shape in a parametric way. Among the places of their possible application are the generation of circles or volumetric objects of rotation, when calculating the Fourier transform, procedural generation of waves on the water plane, generators for a software sound synthesizer, and so on. In all these cases, sin()



and cos()



are called inside the loop, as here:

 for(int n=0; n < num; n++) { const float t = 2.0f*PI*(float)n/(float)num; const float s = sinf(t); const float c = cosf(t); // do something with s and c ... }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

We begin to rewrite the cycle incrementally (see the code below), so it is easier for us to imagine that at iteration n



this cycle with phase t



, the next iteration, n+1



, will consider sin()



and cos()



for t+f



. In other words, sin(t)



and cos(t)



are counted for us and we need to count sin(t+f)



and cos(t+f)



:

 const float f = 2.0f*PI/(float)num; const float t = 0.0f; for( int n=0; n < num; n++ ) { const float s = sinf(t); const float c = cosf(t); // do something with s and c ... t += f; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

It doesn't matter how we got t



and what its range of values ​​is (in the example above - $$$[0; 2 \ pi]$ ) The only thing that bothers us is that there is a loop that constantly calls sin()



and cos()



with a parameter that increases in constant steps (in this case, $$$\ frac {2 \ pi} {\ text {num}}$ ) This article is about how to optimize this code for speed in such a way that the same calculations can be performed without using the sin()



or cos()



functions (in the inner loop), or even the faster sincos()



function.

But if you look at the first formula in the article, we will see that if $$$f = \ alpha$ and $$$t = \ beta$ we can rewrite it as

 sin(t+f) = sin(f)*cos(t) + cos(f)*sin(t) cos(t+f) = cos(f)*cos(t) - sin(f)*sin(t)
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

or in other words

 new_s = sin(f)*old_c + cos(f)*old_s new_c = cos(f)*old_c - sin(f)*old_s
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

Since f



is a constant, sin(f)



and cos(f)



also. Call them a



and b



respectively:

 new_s = b*old_c + a*old_s new_c = a*old_c - b*old_s
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

This expression can be directly used in our code, and then we get a loop in which expensive (and indeed no) trigonometric functions are called!

 const float f = 2.0f*PI/(float)num; const float a = cosf(f); const float b = sinf(f); float s = 0.0f; float c = 1.0f; for( int n=0; n < num; n++ ) { // do something with s and c ... const float ns = b*c + a*s; const float nc = a*c - b*s; c = nc; s = ns; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

Interpretation

To date, we have blindly played with math, not understanding what is really happening. Let's rewrite the inner loop like this:

$$

$$s_ {n + 1} = s_na + c_nb \\ c_ {n + 1} = c_na - s_nb$$

Some of you may have noticed that this is a formula for rotating an object in two-dimensional space. If you still do not understand this, perhaps the matrix form will help you.

$$

$$\ left (\ begin {matrix} s_ {n + 1} \\ c_ {n + 1} \ end {matrix} \ right) = \ left (\ begin {matrix} a & b \\ -b & a \ end {matrix} \ right) \ cdot \ left (\ begin {matrix} s_ {n} \\ c_ {n} \ end {matrix} \ right)$$

In fact, sin(t)



and cos(t)



can be grouped into a vector of length 1 and drawn as a point on the screen. Call this vector x



. Then, $$$x = \ {\ sin \ beta, \ cos \ beta \}$ . So the vector form of expression is $$$x_ {n + 1} = Rx_n$ where $$$R = \ left (\ begin {matrix} a & b \\ - b & a \ end {matrix} \ right)$ . We see that our loop performs a small two-dimensional rotation each iteration so that x



rotates in a circle during the execution of the cycle. Each iteration x



rotates by $$$\ frac {2 \ pi} {\ text {num}}$ radian.

So basically

$$

$$\ sin (\ alpha + \ beta) = \ sin \ alpha \ times \ cos \ beta + \ cos \ alpha \ times \ sin \ beta \\ \ cos (\ alpha + \ beta) = \ cos \ alpha \ times \ cos \ beta - \ sin \ alpha \ times \ sin \ beta$$

this is the point motion formula $$$x = \ {\ cos \ alpha, \ sin \ alpha \}$ circumference in increments of $$$\ beta$ radian. To do this, we will build one of two rotation axes, for example, $$$u = \ {\ cos \ beta, \ sin \ beta \}$ . The first component of the rotation is the projection. $$$x$ on $$$u$ . Because $$$x$ and $$$u$ normalized (have a length of 1), the projection is their scalar product. Consequently, $$$s = x \ cdot u = \ sin \ alpha \ cdot \ cos \ beta + \ cos \ alpha \ cdot \ sin \ beta$ , and of course the second component is the anti-projection, which can be found by projecting onto the perpendicular axis, $$$v$ . We can create this vector by expanding the coordinates $$$u$ and change the sign to the opposite at the first coordinate: $$$c = x \ cdot v = \ cos \ alpha \ cdot \ cos \ beta + \ sin \ alpha \ cdot \ sin \ beta$

Notes

Usually you should be able to perform these tiny rotations over and over again. Indeed, $$$| R | = \ left | \ begin {matrix} a & b \\ - b & a \ end {matrix} \ right | = a ^ 2 + b ^ 2 = \ sin ^ 2 \ alpha + \ cos ^ 2 \ alpha = 1$ , which means that the matrix $$$R$ does not increase or decrease the space to which it is applied, which means that $$$x$ will move in a perfect circle. However, because computers are not accurate, $$$x$ will move in a spiral and eventually coincides with the center of the circle of rotation. I did not have such problems, but I think that they can occur with very large num



, i.e. small turning angles.

Example

In Kindercrasher , the 4096-byte demo from 2008 (screenshot on KDPV), a group of spheres pulsates to the music. For this, I counted the Fourier transform of sound. There are algorithms that do this in real time, for example, FFT . However, my code should fit in a few kilobytes, and I decided to go the other way. Instead of implementing FFT, I wrote DFT by its simple definition. You can check it out on wikipedia.

$$

$$X_k = \ sum_ {n = 0} ^ {N-1} {x_ne ^ {- \ frac {2 \ pi i} {N} kn}} \ quad k = 0,1, \ ldots, N-1$$

My function also takes a 16-bit stereo audio buffer, x



, and returns the first 128 frequencies of the sound spectrum of sound y



. See how the inner loop is organized, the one that runs 4096 times: not a single call to the sin()



or cos()



functions, although in other implementations these calls will be.

 #include <math.h> void iqDFT12( float *y, const short *x ) { for( int i=0; i<128; i++ ) { const float wi = (float)i*(2.0f*3.1415927f/4096.0f); const float sii = sinf( wi ); const float coi = cosf( wi ); float co = 1.0f; float si = 0.0f; float acco = 0.0f; float acsi = 0.0f; for( int j=0; j<4096; j++ ) { const float f = (float)(x[2*j+0]+x[2*j+1]); const float oco = co; acco += co*f; co = co*coi - si*sii; acsi += si*f; si = si*coi + oco*sii; } y[i] = sqrtf(acco*acco+acsi*acsi)*(1.0f/32767.0f); } }