c++
Reinventing the wheel : write your own fast sine…
A couple of days ago I tried to self study the mathematical explanation of pi (π) and ended up with very interesting results and ideas 🙂
First of all lets answer the classic question; What is pi?
Pi or π is a mathematical constant whose value is the ratio of any circle‘s circumference to its diameter in Euclidean space… (Wiki)
From its definition it is pretty easy to find it depending on another mathematical function: sinus.
I can hear that two questions rising up in minds rapidly;
1) Why would I need that? -I don’t know the answer, just curiosity 🙂
2) How? – Answer is in the rest of the post… continue… cmon, go on!
Since it is hard to write mathematical formulas and drawing shapes on our stupid HCI (human computer interaction; keyboard, mouse, etc.) slavery, I prefered to write them on a notebook and took their shots for the ease of understanding and publishing.
Relation Between Pi and Sinus Function
At first galance it is necessery to know that a circle can be expressed by infinite number of triangles,
so let’s start with the most basic equilateral; square.
The first half of the image above shows how to find the area of the square by using its half-diagonal which is the first step to do inductive reasoning.
And the second half is a generalisation of the formula for an equilateral with ‘n’ sides/corners.
As much as we increase the number of sides of the equilateral as much as it approximates to a circle,
Let’s think that we have an equilateral with infinite sides, then it turns into a circle which means in the result formula if we give a very big value to ‘n’ and 1 to ‘r’, then the result approximates to pi.
instead of depending on the number of sides in the equilateral, we want to be dependent to the angle alpha, to do that we simply replace n with 360/alpha.
the conclusion of this part is awesome; we can find the value of sinus in the degree range of zero and ten, with an acceptable error (+0.001) by only one multiplication.
Some Trigonometry
Althoug it is possible to achieve relatively accurate results with the formula shown above for angles between 0 and 10, as the angle gets wider as it loses its accuricy. Therefore we have to use the magic formula for angles less than 10 but how?!
The answer comes from the trigonometric sine addition formula;
sin(a+b) = sin(a) cos(b) + sin(b) cos(a)
If we can keep the ‘b’ less than 10 then we will be able to use our formula in order to find the sine with a couple of aritchmetic operations.
Let’s say we are asked the sine value for 71.654, then;
a = 70
b = 1.654
and,
sin(71.654) = sin(70 + 1.654) = sin(70) cos(1.654) + sin(1.654) cos (70)
In this formula we are able to use the fast calculation for the sin(1.654) part and for the rest unfortunately we need to have sine and cosine tables. The good thing is we only need the multiply of tens for sine and natural number angles between 0 and 10 for cosine. So lets start writing our own fast sine function.
CODING
//precision types
#ifdef PRECISE
#define PRECISION_TYPE double
#else
#define PRECISION_TYPE float
#endif
PRECISION_TYPE hollyConstant = 0.017453292519943295769236907684886;
//First of all sine and cosine tables
PRECISION_TYPE sinTable[] = {
0.0, //sin(0)
0.17364817766693034885171662676931 , //sin(10)
0.34202014332566873304409961468226 , //sin(20)
0.5 , //sin(30)
0.64278760968653932632264340990726 , //sin(40)
0.76604444311897803520239265055542 , //sin(50)
0.86602540378443864676372317075294 , //sin(60)
0.93969262078590838405410927732473 , //sin(70)
0.98480775301220805936674302458952 , //sin(80)
1.0 //sin(90)
};
PRECISION_TYPE cosTable[] = {
1.0 , //cos(0)
0.99984769515639123915701155881391 , //cos(1)
0.99939082701909573000624344004393 , //cos(2)
0.99862953475457387378449205843944 , //cos(3)
0.99756405025982424761316268064426 , //cos(4)
0.99619469809174553229501040247389 , //cos(5)
0.99452189536827333692269194498057 , //cos(6)
0.99254615164132203498006158933058 , //cos(7)
0.99026806874157031508377486734485 , //cos(8)
0.98768834059513772619004024769344 //cos(9)
};
// sin (a+b) = sin(a)*cos(b) + sin(b)*cos(a)
// a = 10*m where m is a natural number and 0<= m <= 90
// i.e. lets a+b = 18.22
// then a = 10, b = 8.22
PRECISION_TYPE myFastSin ( PRECISION_TYPE angle )
{
int a = angle * 0.1f;
PRECISION_TYPE b = angle - 10 * a;
return sinTable[a] * cosTable[int(b)] + b * hollyConstant * sinTable[9-a];
}
RESULTS and ERROR rate
The results are really astonishing 😛 I wasn’t really expecting such a result but it just did happen 😀
the first result is the result of our fast sine function and the second one is the result of <cmath> ‘s result
0. angle = 0.803589, expected = 0.0140248
OK : FASTER !! time dif. = 4294967233
result = 0.0140253, time spent = 29, error = -4.61005e-007, error % = -0.00328707
result = 0.0140248, time spent = 92, error = 0, error % = 01. angle = 1.50818, expected = 0.0263196
OK : FASTER !! time dif. = 4294967262
result = 0.0263227, time spent = 31, error = -3.03984e-006, error % = -0.0115497
result = 0.0263196, time spent = 65, error = 0, error % = 02. angle = 2.21377, expected = 0.0386279
OK : FASTER !! time dif. = 4294967261
result = 0.0386375, time spent = 31, error = -9.61125e-006, error % = -0.0248816
result = 0.0386279, time spent = 66, error = 0, error % = 03. angle = 2.92035, expected = 0.0509477
OK : FASTER !! time dif. = 4294967259
result = 0.0509698, time spent = 29, error = -2.20649e-005, error % = -0.0433089
result = 0.0509477, time spent = 66, error = 0, error % = 04. angle = 3.62794, expected = 0.0632772
OK : FASTER !! time dif. = 4294967261
result = 0.0633195, time spent = 30, error = -4.23118e-005, error % = -0.0668674
result = 0.0632772, time spent = 65, error = 0, error % = 05. angle = 4.33653, expected = 0.0756145
OK : FASTER !! time dif. = 4294967259
result = 0.0756867, time spent = 30, error = -7.22408e-005, error % = -0.0955383result = 0.0756145, time spent = 67, error = 0, error % = 0
….
94. angle = 71.4059, expected = 0.947801
OK : FASTER !! time dif. = 4294967260
result = 0.947942, time spent = 30, error = -0.000140667, error % = -0.0148414
result = 0.947801, time spent = 66, error = 0, error % = 095. angle = 72.2045, expected = 0.952153
OK : FASTER !! time dif. = 4294967263
result = 0.95228, time spent = 30, error = -0.000126302, error % = -0.0132649
result = 0.952153, time spent = 63, error = 0, error % = 096. angle = 73.0041, expected = 0.956326
OK : FASTER !! time dif. = 4294967262
result = 0.956337, time spent = 29, error = -1.17421e-005, error % = -0.00122784
result = 0.956326, time spent = 63, error = 0, error % = 097. angle = 73.8047, expected = 0.960316
OK : FASTER !! time dif. = 4294967264
result = 0.961116, time spent = 31, error = -0.000799954, error % = -0.0833011
result = 0.960316, time spent = 63, error = 0, error % = 098. angle = 74.6063, expected = 0.964124
OK : FASTER !! time dif. = 4294967261
result = 0.9649, time spent = 29, error = -0.000775695, error % = -0.0804559
result = 0.964124, time spent = 64, error = 0, error % = 099 times faster
0 times more precise
%218.04 faster
total error % = 9.46265, maximum error = 0.00242305, minimum error = 4.47035e-008, average error % =0.0955823
so the average error is %0.1 and its more than 2 times faster. On the other hand the memory used for the sine and cosine tables also should be taken into account but I don’t think its a big deal anyway.