strawberry
a sin/cos trick



You probably remember this identity from school:



When you are 14 and are first introduced to it you probably cannot do much about it but memorize it by heart. Too bad, cause in fact you don't need to do that - this formula has a very concise visual interpretation. In fact, when I need to write the formula above I just make a drawing, analize it, and write the formula as I interpret the image. This image will be obvious by the time we reach half this tutorial. But for now, to keep the things fun and delay the euroeka moment, let's first think why we would care about this formula at all.


Intro


The sin and cos trigonometric functions are probably the functions that most frequently appear in computer graphics, for they are the basics to describe any circular shape in a parametric way. Among their uses, we have the generation of circles or 3D revolution objects, when computation of a fourier transform, procedural waves for a water plane, generators for a software sound synthetizer, etc etc. In all these cases, sin, cos or both are repeatedly called inside a loop, similar to this:

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

}
Kindercrasher (2008), using the trick explained here for eficiently extracing sound frequency information
For the point of this articles the exact way "t" is computed doesn't matter, neither does the range of values it takes (0 to 2PI in this example above). The only thing we are concerned with is that there is a loop repeatedly calling sin and cos with a parameter that increments in constant steps (2·PI/num, in this example). This article is about how to optimize this code, for speed, such that the same computations can be performed without using the sin or cos functions at all (in the inner loop), not even the faster sincos() function.



The transformation


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;
}
incremental version of the loop 		We start by rewritting our loop in an incremental way (see code on the left), so that it easier to see that given an iteration n of the loop which has phase t, the next iteration n+1 is gonna compute the sin and cos of t+f. In other words, given that we just computed sin(t) and cos(t), we now have to compute sin(t+f) and cos(t+f). But looking to the formula above, we see that if we set f to be alpha, t to be beta, we can rewrite this 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, so are cos(f) and sin(f). Let's call them a and b:

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

This expression can be directly used in out code, resulting in a loop that hasn't any expensive trigonometric function at all!

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


The interpretation


So far we have blindly played with a math indentity without really seeing what we were doing. Let's first rewrite the inner loop this way:



Some of you might have already noticed that this is nothing but the formula for a 2D rotation. If you didn't recognize it yet, perhaps the matricial form of it might help:



Indeed, sin(t) and cos(t) can be grouped as a vector of length one (and ploted as a dot in the screen). Let's call it x:



Then, the vectorial form of the expression above looks like



with R being, of course, the {a, b, -b, a} 2x2 matrix. So, we see that out loop is performing a little 2D rotation in each iteration, such that x rotates in a circle during the execution of the loop. Every iteration x rotates by an angle of 2·PI/num radians.
So, basically



is a 2D rotation formula of a point x = {cos(alpha), sin(alpha)} in the circle by an angle of beta radians. To do so, we first contruct one of the two axes of the rotation, ie, u = { cos(beta), sin(beta) }. The first component of the rotation is gonna be the projection of x into u. Cause both x and u have length 1, the projection is just their dot product. Therefore:



and of course, the second component is it's antiprojection, which we can do by projecting in a perpendicular axis, v. We can construct this vector by reversing the coordinates of u and negating one of them:





Some notes


Normally you should be able to perform this tiny rotations over and over again. Indeed,



what means that R is not shrinking nor expanding the space as it is applied, or in other words, it means that x will move in a perfect cicle. However due to numerical precission issues, x move in a spiral and fall into the origin. I never had this issue in my applications, but I'm guessing that for very big values of num, ie, little rotation angles, one can get some problems.



An example


In Kindercrasher, a realtime 4 kilobytes demo from 2008, a group of spheres do pulsate to the music. For that I computed the Fourier transform of the sound. Fast algorithms exist to do so in realtime, like the FFT. However, given that my code had to fit in no more that few hundred bytes, I decided to do it in a different way. Instead of implementing the FFT, I implemented the Discrete Fourier Transform from it's very basic definition. You can check it on the wikipedia.



My function takes an stereo 16bit sound buffer, x, and returns the first 128 frequencies in the sound spectrum of the sound in y. Note how the inner loop, that one iterating 4096 times, has no sin or cos calls, as the naive implementation would. 
#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);
    }
}




Back


iņigo quilez 2010