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 and 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 recursively applying 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 rotation 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 determinant as:

The new transformation formulas for are:

Which are more complex than the original formulas (though they work now).
It’s also possible to further improve them noticing that appears in computation.
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 and 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’..)