## 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 applying$R'$ 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 of$s$, 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 within$R'$, 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 write$R''$ 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 that$u'_x$ appears in$u'_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 devote$n$ 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); } 
(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:

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)