Inigo Quilez   ::     ::  
You probably remember this identity from school:

sin(α+β)=sin α ⋅ cos β + cos α ⋅ sin β
cos(α+β)=cos α ⋅ cos β - sin α ⋅ sin β

When you are kid and are first introduced to it your first feeling is the pain of having to memorize it. That's too bad, because in fact you don't need to memorize anthing - this formula arises naturally when you draw the rotation of a triangle on a piece of paper. In fact, that's exactly what I do when I need to write this formula down. This interpretation will be obvious by the time we are half way down this tutorial. But for now, to keep the things fun and delay the eureka moment, let's first think why we would care about this formula at all.


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 synthesizer, 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 ... }

We start by rewriting 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):

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

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.

But looking to the first formula in this article, the one we memorized a kids, 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 our 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 identity without really seeing what we were doing. Let's first rewrite the inner loop this way:

sn+1 = sn ⋅ a + cn⋅b
cn+1 = cn ⋅ a - sn⋅b

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 plotted as a dot in the screen). Let's call it x:

x = { sin β, cos β }

Then, the vectorial form of the expression above looks like

xn+1 = R ⋅ xn

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π/num radians.
So, basically

sin(α+β)=sin α ⋅ cos β + cos α ⋅ sin β
cos(α+β)=cos α ⋅ cos β - sin α ⋅ sin β

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

s = <x,u> = sin α ⋅ cos β + cos α⋅ sin β

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:

c = <x,v> = cos α ⋅ cos β - sinα⋅ sin β

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 circle. However due to numerical precision 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); } }