When I was writing the Nintendo DS 4k intro, I obviously needed sin and cos functions to set up rotation matrices. Since using the C standard library for this task would have exceeded the 4kb limit easily, I was looking for an alternative that does not use much memory.
I was searching for an approximation method using taylor series and google made me visit these sites:
- polygonal labs – Fast and accurate sine/cosine approximation
- devmaster.net – Fast and accurate sine/cosine
The only task left for me was porting it from float to fixed point. I didn’t really care about precision or performance much, something that reminds to a sine-wave and isn’t horribly slow would suffice.
Suprisingly, the fixed-point sin approximation is quite precise. At least more precise than I though it would be. Here is a screenshot of the original sinf routine (green) and the approximated one (red) using the source code below:
Well, there is an issue with my code. When the incoming radians parameter is getting bigger (1000 and above), precision is getting worse, but this is an issue that can be fixed if you’re up to.
The routine from polygon labs can’t really handle radians greater than 2*PI, which I tried to fix by normalizing it to an angle between -1..+1 and then just use a bit-AND to discard all numbers outside this range.
I wanted to keep the radians format for the incoming parameter, because it’s how the standard sinf routine works. Using another range, like 0..4096 to represent a full rotation would probably solve that precision problem with large radians because this would remove the need for the first two code lines, but I want 2*PI.
Here is the source code, to be compiled with libnds and devkitARM:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
// Gets sin of specified radians. // Input and output format is fixed point sign.19.12 int32 f32sin(int32 radians) { // Offset to slightly restore precision when radians is large. const int32 offset= f32mul(radians, floattof32(0.00147f)); // Divide incoming radians by 2*PI, so it is a "normalized" // angle between -1 and +1, where 1 equals 2*PI (360 degree). int32 value = f32mul(radians + offset, floattof32(1.0f / (2*M_PI))); // Now that we have our normalized angle, // we just discard all numbers which are not in the range -1..+1. // Since this is fixed point, a simple AND operation can be used. value &= floattof32(1)-1; // Always wrap normalized angle to -0.5(-PI)..+0.5(+PI) if(value > floattof32(0.5f)) value -= floattof32(1); // subtract 2*PI in normalized form else if(value < floattof32(-0.5f)) value += floattof32(1); // add 2*PI in normalized form // Convert normalized angle (-1..1) back to radians (-2*PI..2*PI) value = f32mul(value, floattof32(2*M_PI)); // Approximate sin value const int32 B = floattof32(1.2732395447351f); const int32 C = floattof32(-0.405284734569f); return f32mul(B, value) + f32mul(f32mul(C, value), f32abs(value)); } inline int32 f32abs(int32 a) { return a >= 0 ? a : -a; } // multiplies two fixed point numbers. // since the result is shifted rather than the two // operands, make sure a and b are small. inline int32 f32mul(int32 a, int32 b) { return (a*b)>>12; } |
Let me know when you fixed the precision problem with large radian values 🙂
Edit: Jasper “cearn” Vijn apparently likes challenges or in his own words:
So here I am, looking forward to a nice quiet weekend; hang back, watch some telly and maybe read a bit – but NNnnneeeEEEEEUUUuuuuuuuu!! Someone had to write an interesting article about sine approximation. With a challenge at the end. And using an inefficient kind of approximation. And so now, instead of just relaxing, I have to spend my entire weekend and most of the week figuring out a better way of doing it. I hate it when this happens >_<.
Don’t miss to read his article Another fast fixed-point sine approximation.
For completeness, I post a fixed version of the method above. I highly encourage you to visit Jasper’s site and get one of his better sine approximation functions though.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
// Gets sin of specified radians. // Input and output format is fixed point sign.19.12 int32 f32sin(int32 radians) { const int32 PI2_PRECISION_BITS = 20; const int32 PI2_RECIPROCAL = (1.0f / (2*M_PI)) * (1<<PI2_PRECISION_BITS); // Divide incoming radians by 2*PI, so it is a "normalized" // angle between -1 and +1, where 1 equals 2*PI (360 degree). int32 value = (radians * PI2_RECIPROCAL) >> PI2_PRECISION_BITS; // Now that we have our normalized angle, // we just discard all numbers which are not in the range -1..+1. // Since this is fixed point, a simple AND operation can be used. value &= floattof32(1)-1; // Always wrap normalized angle to -0.5(-PI)..+0.5(+PI) if(value > floattof32(0.5f)) value -= floattof32(1); // subtract 2*PI in normalized form else if(value < floattof32(-0.5f)) value += floattof32(1); // add 2*PI in normalized form // Convert normalized angle (-1..1) back to radians (-2*PI..2*PI) value = f32mul(value, floattof32(2*M_PI)); // Approximate sin value const int32 B = floattof32(1.2732395447351f); const int32 C = floattof32(-0.405284734569f); return f32mul(B, value) + f32mul(f32mul(C, value), f32abs(value)); } |