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
within a certain interval, using infinitesimal rotations.We start to work in two dimensions with a clockwise rotation matrix:
For very small rotation angles
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:
Now we can generate our values on an interval
on a unit vector
Supposing that we don’t have at hand any code or library that can compute
for an arbitrary value of
, we start with
Unfortunately this approach is not very good: if we plot
, iteration after iteration, we observe that after a full
is not a unit vector anymore.
Additional rotations will typically plot a diverging spiral as gets bigger and bigger.
Numerical issues aside the problem lies within
, which is not a rotation matrix anymore; in fact its determinant is:
Areas are not conserved, no matter how accurate our computation is. Since we want
to represent a proper rotation matrix we can try to fix it in order to have a new rotation matrix
with an unitary determinant.
We arbitrarily write
The new transformation formulas for
Which are more complex than the original formulas (though they work now).
It’s also possible to further improve them noticing that
The final transformation as:
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
possible is equal to
if we devote
bits to the fractional part of our fixed point numbers, then just
integer additions and subtractions are needed to generate all the values of
within the interval
, 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’..)