A simple way to pre-compute sin/cos tables

This is a kind of an old school trick, rarely useful nowadays, but it’s nonetheless interesting due to its simplicity and elegancy.
The main idea is about filling a table with pre-computed values of\sin a and\cos a within a certain interval, using infinitesimal rotations.We start to work in two dimensions with a clockwise rotation matrix:

R = \begin{pmatrix}\cos\alpha\ && -\sin\alpha\  \\ \sin\alpha\ && \cos\alpha\ \end{pmatrix}

For very small rotation angles\alpha\ \to 0 thus we can expand our rotation matrix elements with a Taylor series around zero.
Ignoring all the non linear terms we get something like this:

R' = \begin{pmatrix}1  && -\alpha\  \\ \alpha\  && 1 \end{pmatrix}

Now we can generate our values on an interval\left [ s, t \right ) recursively applyingR' on a unit vector\hat{u} = ( \cos s, \sin s ) .

\hat{u'} = \hat{u}R'

=( u_x, u_y ) \begin{pmatrix}1  && -\alpha\  \\ \alpha\ &&  1 \end{pmatrix}

= \begin{cases} u'_x = u_x +  u_y\alpha\ \\ u'_y = -u_x\alpha\ + u_y \end{cases}

Supposing that we don’t have at hand any code or library that can compute\hat{u} for an arbitrary value ofs , we start with s = 0 \to \hat{u} = ( 1, 0 ).

Unfortunately this approach is not very good: if we plot\hat{u} , iteration after iteration, we observe that after a full 2\pi rotation \hat{u} is not a unit vector anymore.
Additional rotations will typically plot a diverging spiral as \big\| \overrightarrow{u} \big\| gets bigger and bigger.
Numerical issues aside the problem lies withinR' , which is not a rotation matrix anymore; in fact its determinant is:

\begin{vmatrix}1  && -\alpha\  \\ \alpha\  && 1 \end{vmatrix} = 1 + \alpha\ ^2 \ne 1

Areas are not conserved, no matter how accurate our computation is. Since we want R' to represent a proper rotation matrix we can try to fix it in order to have a new rotation matrix R''  with an unitary determinant.
We arbitrarily writeR''  determinant as:

\begin{vmatrix}1  && -\alpha\  \\ \alpha\ && 1 - \alpha\ ^2 \end{vmatrix} = 1

The new transformation formulas for\hat{u} are:

\begin{cases} u'_x = u_x +  u_y\alpha\ \\ u'_y = -u_x\alpha\ + u_y - u_y\alpha\ ^2 \end{cases}

Which are more complex than the original formulas (though they work now).
It’s also possible to further improve them noticing thatu'_x appears inu'_y computation.
The final transformation as:

\begin{cases} u'_x = u_x +  u_y\alpha\ \\ u'_y = -u'_x\alpha\ + u_y \end{cases}

These final formulas are unexpectedly simple and also much more accurate than the previous ones, also using fixed point math we get that the smallest\alpha\ possible is equal to 1 .
if we devoten bits to the fractional part of our fixed point numbers, then just 2^{n + \pi} integer additions and subtractions are needed to generate all the values of\sin\alpha\ and\cos\alpha\ within the interval\left [ 0, 2\pi \right ) , unless we want to perform less iterations and then use symmetry properties to generate the missing values.
Here some example code (not really meant to be fast or anything):

const int32_t fracBits = 12;
const int32_t iterationCount = (int16_t)(2.0f * 3.4159f * (1<<fracBits));
static int32_t cosTable[iterationCount];
static int32_t sinTable[iterationCount];

int32_t tempCos = 1<<(fracBits<<1);
int32_t tempSin = 0;

for (uint32_t i = 1; i < iterationCount; i++)
    cosTable[i] = tempCos>>fracBits;
    sinTable[i] = tempSin>>fracBits;
    tempCos = tempCos + (tempSin>>fracBits);
    tempSin = tempSin - (tempCos>>fracBits);

That’s it for now, any comments and/or advices are welcome!
(and pls don’t be to picky about the funny text and formulas formatting, I’m still learning how to write on this ‘thing’..)


5 Responses to “A simple way to pre-compute sin/cos tables”

  1. Arun Demeure Says:

    Hey Marco, good luck on the blog, hopefully there will be plenty of cool stuff here soon!

    That is a very elegant technique indeed – wouldn’t have imagined you could build a sin/cos table that efficiently. It’s obviously overkill for PCs and modern consoles but I wonder if it could be a good optimization on handhelds. A Nintendo DS isn’t exactly super-fast…

    On a modern PC/console, my preference for sin/cos tends to be the combination of these two techniques, the first is from Nick (also on the Beyond3D forums) and the second is from Christer Ericson, who I notice already commented in the other post! 🙂


  2. marco Says:

    Ciao Arun, thanks for your words, I’ll try to do everything I can do to make this blog interesting to read.
    Good idea about the DS, I didn’t think about it, though I have no idea what kind of floating point capabilities that little marvel has..

  3. Arun Demeure Says:

    Googling around quickly tells me it has no FPU at all… So AFAICT, this technique could indeed be quite nifty there.

  4. christerericson Says:

    Marco, the code looks like you transscribed it wrong, because if you substitute the first line of the loop body for the “cosTable[i+1]” in the second line, you get “sinTable[i+1] = -cosTable[i]” which can’t be right.

    Anyway, I was going to say that this is somewhat similar to Goertzel’s algorithm: http://www.dattalo.com/technical/theory/sinewave.html

    You can also get similar code from the “recursive” identities:

    sin(a+b) = sin(a)cos(b) + cos(a)sin(b)
    cos(a+b) = cos(a)cos(b) – sin(a)sin(b)

    The latter approach is also described in the link above, in section 4.

  5. marco Says:

    Hello Christer,

    you’re obviously right, I forgot to renormalize my fixed point quantities.
    The code is fixed now, it’s still cheap but there are 4 additional shifts and hopefully no additional bugs!
    I also substituted the word ‘versor’ with unit vector, according Jim Tilander only in Italy we use that word (ciao Jim!) 🙂

    I’ve never heard about Goertzel’s algorithm before, but thanks for yuor link, it’s full of useful tricks (the page about logarithm approximation is quite interesting)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: