math: Sine and cosine accuracy

[mansourmoufid]

Hi,

When I try to test the sine and cosine implementations with inputs on [0, 2pi] rather than [-1, 1], they are only accurate to two or three decimal places.

I know people put in the effort and I don't want to be rude and mess with their code, but I would like to suggest two formulas which are more accurate and probably just as fast. They are polynomial approximations accurate to nine decimal places. The formulas are here: http://eliteraspberries.github.io/pal/formulas/

Polynomial approximations involve mainly floating-point multiplication and addition, which is good for modern hardware with a multiply-accumulate instruction.

[DanXS]

I can confirm that the accuracy of the sin function slips as the value of theta gets larger, say greater than pi. It could be related to symmetry. I ran a test using just a c version on my computer as I don't have a parallela. I've also provided an alternate approach which is slightly more memory hungry as it uses a lookup table and bilinear interpolation. I ran a test for a smooth range from [0 to 2pi], the results show the operating systems sin function in comparison to the iterative method implemented here and my approach of using interpolation and some pre-calculated values. In this example I pre-calculated 360 values to match the standard number of degrees although the table could be larger or smaller if required. The other values are calculated via bilinear interpolation. Note that there are divide operations in the code but only for constants which should be removed by the compiler so wouldn't need to be calculated on the Parallela. It is really a matter of whether it is worth using a little extra memory and getting accurate results or whether the iterative method can be improved by taking advantage of symmetries or perhaps more iterations?

I've attached the test code - please build and test for yourselves:

`

#include <stdio.h>
#include <math.h>

const double sinTable[] = {
0.0, 0.01745240643728351, 0.03489949670250097, 0.05233595624294383, 0.0697564737441253,
0.08715574274765817, 0.10452846326765346, 0.12186934340514748, 0.13917310096006544, 0.15643446504023087,
0.17364817766693033, 0.1908089953765448, 0.20791169081775931, 0.224951054343865, 0.24192189559966773,
0.25881904510252074, 0.27563735581699916, 0.29237170472273677, 0.3090169943749474, 0.32556815445715664,
0.3420201433256687, 0.35836794954530027, 0.374606593415912, 0.3907311284892737, 0.40673664307580015,
0.42261826174069944, 0.4383711467890774, 0.45399049973954675, 0.4694715627858908, 0.48480962024633706,
0.49999999999999994, 0.5150380749100542, 0.5299192642332049, 0.5446390350150271, 0.5591929034707469,
0.573576436351046, 0.5877852522924731, 0.6018150231520483, 0.6156614753256582, 0.6293203910498374,
0.6427876096865393, 0.6560590289905072, 0.6691306063588582, 0.6819983600624985, 0.6946583704589973,
0.7071067811865475, 0.7193398003386511, 0.7313537016191705, 0.7431448254773941, 0.7547095802227719,
0.766044443118978, 0.7771459614569708, 0.7880107536067219, 0.7986355100472928, 0.8090169943749473,
0.8191520442889917, 0.8290375725550416, 0.8386705679454239, 0.848048096156426, 0.8571673007021121,
0.8660254037844386, 0.8746197071393957, 0.8829475928589269, 0.8910065241883678, 0.898794046299167,
0.9063077870366499, 0.9135454576426009, 0.9205048534524403, 0.9271838545667874, 0.9335804264972017,
0.9396926207859083, 0.9455185755993167, 0.9510565162951535, 0.9563047559630354, 0.9612616959383189,
0.9659258262890683, 0.9702957262759965, 0.9743700647852352, 0.9781476007338056, 0.981627183447664,
0.984807753012208, 0.9876883405951378, 0.9902680687415703, 0.992546151641322, 0.9945218953682733,
0.9961946980917455, 0.9975640502598242, 0.9986295347545738, 0.9993908270190958, 0.9998476951563913,
1.0, 0.9998476951563913, 0.9993908270190958, 0.9986295347545738, 0.9975640502598242, 0.9961946980917455,
0.9945218953682734, 0.9925461516413221, 0.9902680687415704, 0.9876883405951377, 0.984807753012208,
0.981627183447664, 0.9781476007338057, 0.9743700647852352, 0.9702957262759965, 0.9659258262890683,
0.9612616959383189, 0.9563047559630355, 0.9510565162951536, 0.9455185755993168, 0.9396926207859084,
0.9335804264972017, 0.9271838545667874, 0.9205048534524404, 0.913545457642601, 0.90630778703665,
0.8987940462991669, 0.8910065241883679, 0.8829475928589271, 0.8746197071393959, 0.8660254037844388,
0.8571673007021123, 0.8480480961564261, 0.838670567945424, 0.8290375725550418, 0.819152044288992,
0.8090169943749475, 0.7986355100472928, 0.788010753606722, 0.7771459614569711, 0.7660444431189781,
0.7547095802227719, 0.7431448254773942, 0.7313537016191707, 0.7193398003386514, 0.7071067811865476,
0.6946583704589971, 0.6819983600624985, 0.6691306063588583, 0.6560590289905073, 0.6427876096865395,
0.6293203910498377, 0.6156614753256584, 0.6018150231520482, 0.5877852522924732, 0.5735764363510464,
0.5591929034707469, 0.5446390350150269, 0.5299192642332049, 0.5150380749100544, 0.49999999999999994,
0.48480962024633717, 0.4694715627858911, 0.45399049973954686, 0.4383711467890773, 0.4226182617406995,
0.40673664307580043, 0.39073112848927416, 0.37460659341591224, 0.3583679495453002, 0.3420201433256689,
0.325568154457157, 0.3090169943749475, 0.292371704722737, 0.27563735581699966, 0.258819045102521,
0.24192189559966773, 0.22495105434386475, 0.2079116908177593, 0.19080899537654497, 0.17364817766693025,
0.15643446504023098, 0.13917310096006574, 0.12186934340514755, 0.10452846326765373, 0.08715574274765864,
0.06975647374412552, 0.05233595624294381, 0.0348994967025007, 0.017452406437283435, 1.2246467991473532e-16,
-0.017452406437283192, -0.03489949670250089, -0.05233595624294356, -0.06975647374412483, -0.08715574274765794,
-0.10452846326765305, -0.12186934340514774, -0.13917310096006552, -0.15643446504023076, -0.17364817766693047,
-0.19080899537654475, -0.20791169081775907, -0.22495105434386498, -0.2419218955996675, -0.25881904510252035,
-0.275637355816999, -0.2923717047227364, -0.30901699437494773, -0.3255681544571568, -0.34202014332566866,
-0.35836794954530043, -0.374606593415912, -0.39073112848927355, -0.4067366430757998, -0.4226182617406993,
-0.43837114678907707, -0.4539904997395463, -0.46947156278589086, -0.48480962024633695, -0.5000000000000001,
-0.5150380749100542, -0.5299192642332048, -0.5446390350150271, -0.5591929034707467, -0.5735764363510458,
-0.587785252292473, -0.601815023152048, -0.6156614753256578, -0.6293203910498376, -0.6427876096865393,
-0.6560590289905075, -0.6691306063588582, -0.6819983600624984, -0.6946583704589974, -0.7071067811865475,
-0.7193398003386509, -0.7313537016191701, -0.743144825477394, -0.7547095802227717, -0.7660444431189779,
-0.7771459614569711, -0.788010753606722, -0.7986355100472928, -0.8090169943749473, -0.8191520442889916,
-0.8290375725550414, -0.838670567945424, -0.8480480961564258, -0.8571673007021121, -0.8660254037844384,
-0.8746197071393959, -0.882947592858927, -0.8910065241883678, -0.8987940462991668, -0.9063077870366497,
-0.913545457642601, -0.9205048534524403, -0.9271838545667873, -0.9335804264972016, -0.9396926207859082,
-0.9455185755993168, -0.9510565162951535, -0.9563047559630353, -0.961261695938319, -0.9659258262890683,
-0.9702957262759965, -0.9743700647852351, -0.9781476007338056, -0.9816271834476639, -0.984807753012208,
-0.9876883405951377, -0.9902680687415704, -0.9925461516413221, -0.9945218953682734, -0.9961946980917455,
-0.9975640502598242, -0.9986295347545738, -0.9993908270190957, -0.9998476951563913, -1.0,
-0.9998476951563913, -0.9993908270190958, -0.9986295347545738, -0.9975640502598243, -0.9961946980917455,
-0.9945218953682734, -0.992546151641322, -0.9902680687415704, -0.9876883405951378, -0.9848077530122081,
-0.9816271834476641, -0.9781476007338058, -0.9743700647852352, -0.9702957262759966, -0.9659258262890682,
-0.9612616959383188, -0.9563047559630354, -0.9510565162951536, -0.945518575599317, -0.9396926207859085,
-0.9335804264972021, -0.9271838545667874, -0.9205048534524405, -0.9135454576426008, -0.9063077870366498,
-0.898794046299167, -0.891006524188368, -0.8829475928589271, -0.8746197071393961, -0.8660254037844386,
-0.8571673007021123, -0.8480480961564262, -0.8386705679454243, -0.8290375725550421, -0.8191520442889918,
-0.8090169943749476, -0.7986355100472932, -0.7880107536067218, -0.7771459614569708, -0.7660444431189781,
-0.7547095802227722, -0.7431448254773947, -0.731353701619171, -0.7193398003386517, -0.7071067811865477,
-0.6946583704589976, -0.6819983600624983, -0.6691306063588581, -0.6560590289905074, -0.6427876096865396,
-0.6293203910498378, -0.6156614753256588, -0.6018150231520483, -0.5877852522924732, -0.5735764363510464,
-0.5591929034707473, -0.544639035015027, -0.5299192642332058, -0.5150380749100545, -0.5000000000000004,
-0.48480962024633684, -0.4694715627858908, -0.45399049973954697, -0.43837114678907696, -0.4226182617407,
-0.40673664307580015, -0.39073112848927466, -0.3746065934159123, -0.35836794954530077, -0.34202014332566855,
-0.32556815445715753, -0.3090169943749476, -0.29237170472273627, -0.2756373558169998, -0.2588190451025207,
-0.24192189559966787, -0.22495105434386534, -0.20791169081775987, -0.19080899537654467, -0.17364817766693127,
-0.15643446504023112, -0.13917310096006588, -0.12186934340514811, -0.10452846326765342, -0.08715574274765832,
-0.06975647374412476, -0.05233595624294437, -0.03489949670250082, -0.01745240643728445};

const double RAD_TO_DEG = 180.0/3.14159265359;
const double DEG_TO_RAD = 1.0/RAD_TO_DEG;

#define wrap(i) (i % 360 < 0) ? (360 + ((i) % 360)) : ((i) % 360)
#define sinLookup(i) (sinTable[wrap(i)])

double interpolated_sin(double theta) {
// convert radians to degrees for look up table
double degrees = RAD_TO_DEG*theta;

// let i be the rounded index
int i = (int)degrees;

// sample from sin lookup table
double p0 = sinLookup(i-1);
double p1 = sinLookup(i);
double p2 = sinLookup(i+1);
double p3 = sinLookup(i+2);

// cubic interpolation coefficient calculations
double a = (-1.0/2.0)*p0 + (3.0/2.0)*p1 - (3.0/2.0)*p2 + (1.0/2.0)*p3;
double b = p0 - (5.0/2.0)*p1 + 2.0*p2 - (1.0/2.0)*p3;
double c = (-1.0/2.0)*p0 + (1.0/2.0)*p2;
double d = p1;

// calculate x, x^2 and x^3
double x = degrees - (double)i;
double x2 = x*x;
double x3 = x*x2;

// return interpolated result
return (a*x3 + b*x2 + c*x + d);

}

double iterative_sin(double theta) {
double theta2 = theta * theta;
double val = 0.0;
val = (1.0) - theta2 * (0.083333333) * (0.076923077)* val;
val = (1.0) - theta2 * (0.1) * (0.090909091) * val;
val = (1.0) - theta2 * (0.125) * (0.111111111) * val;
val = (1.0) - theta2 * (0.166666667) * (0.142857143) * val;
val = (1.0) - theta2 * (0.25) * (0.2) * val;
val = (1.0) - theta2 * (0.5) * (0.333333333) * val;
return theta * val;
}

int main() {

const double steps = 1000;
for (int i = 0; i <= (int)steps; i++) {
    double theta = i/steps*360.0*DEG_TO_RAD;
    printf("sin(%g) = %g\n", theta, sin(theta));
    printf("iterative sin(%g) = %g - error=%g\n", theta, iterative_sin(theta), iterative_sin(theta)-sin(theta));
    printf("interpolated sin(%g) = %g - error=%g\n", theta, interpolated_sin(theta), interpolated_sin(theta)-sin(theta));
}

return 0;

}

`

PS. Is that deal still on where you get a Parallela if you contribute a function - I'm unemployed at the moment!

[DanXS]

I've attached the results and test code here:
sin_test.txt

sin.c.txt