Why does this sin warm the table inaccurately when the radian is large?

Question

I want to create a sin cos lookup table for optimization using an array index from 0 to UCHAR_MAX , so 0 rad is the index 0, pi/2 UCHAR_MAX/4 :

sincos.h

 #include <limits.h> #include <math.h> int sini[UCHAR_MAX]; int cosi[UCHAR_MAX]; #define MAGNIFICATION 256 #define SIN(i) sini[i]/MAGNIFICATION #define COS(i) cosi[i]/MAGNIFICATION void initTable(){ for(int i=0;i<UCHAR_MAX;i++){ sini[i]=sinf(i*2*M_PI/UCHAR_MAX)*MAGNIFICATION; cosi[i]=cosf(i*2*M_PI/UCHAR_MAX)*MAGNIFICATION; } } 

the reason for using UCHAR_MAX as max is that I want to make good use of the unsigned char overflow to simulate a radian that only changes from 0 to 2*pi : for example, if the radian value is 2*pi , the array index becomes UCHAR_MAX because it overflows, it automatically becomes 0, and no mod is required (if I use 0 to 360 as a domain, I may need to calculate index%360 each time). Then I test it with some radians:

 float rad[]={2.0f,4.0f,6.0f,8.0f,10.0f,-2.0f,-4.0f,-6.0f,-8.0f,-10.0f}; 

as below:

 #include "sincos.h" #include <stdio.h> int main(){ initTable(); unsigned char radToIndex; float rad[]={2.0f,4.0f,6.0f,8.0f,10.0f,-2.0f,-4.0f,-6.0f,-8.0f,-10.0f}; int scalar=123; printf("scalar=%d\n",scalar); for(int i=0;i<sizeof(rad)/sizeof(float);i++){ radToIndex=rad[i]*UCHAR_MAX/2/M_PI; printf("%d*sin(%f) : %f , %d\n",scalar,rad[i],scalar*sinf(rad[i]),scalar*SIN(radToIndex)); } return 0; } 

I test the table with 123*sin(radian) , I find that the results begin to go beyond the actual when the radian increases (when the radian is 10 or -10):

 scalar=123 123*sin(2.000000) : 111.843582 , 111 123*sin(4.000000) : -93.086708 , -92 123*sin(6.000000) : -34.368107 , -35 123*sin(8.000000) : 121.691063 , 122 123*sin(10.000000) : -66.914597 , -61 123*sin(-2.000000) : -111.843582 , -112 123*sin(-4.000000) : 93.086708 , 90 123*sin(-6.000000) : 34.368107 , 38 123*sin(-8.000000) : -121.691063 , -122 123*sin(-10.000000) : 66.914597 , 59 

and check other data:

 float rad[]={0.01f,0.1f,1.0f,10.0f,100.0f,1000.0f,-0.01f,-0.1f,-1.0f,-10.0f,-100.0f,-1000.0f}; 

output:

 scalar=123 123*sin(0.010000) : 1.229980 , 0 123*sin(0.100000) : 12.279510 , 12 123*sin(1.000000) : 103.500931 , 102 123*sin(10.000000) : -66.914597 , -61 123*sin(100.000000) : -62.282974 , -97 123*sin(1000.000000) : 101.706184 , -25 123*sin(-0.010000) : -1.229980 , 0 123*sin(-0.100000) : -12.279510 , -8 123*sin(-1.000000) : -103.500931 , -100 123*sin(-10.000000) : 66.914597 , 59 123*sin(-100.000000) : 62.282974 , 98 123*sin(-1000.000000) : -101.706184 , 22 

The error increases when the value increases, so I am sure that the table becomes inaccurate when the radian is large. In sincos.h there is a MAGNIFICATION value for accuracy control, I changed it from 256 to 4096, but it does not seem to improve:

 scalar=123 123*sin(0.010000) : 1.229980 , 0 123*sin(0.100000) : 12.279510 , 12 123*sin(1.000000) : 103.500931 , 102 123*sin(10.000000) : -66.914597 , -62 123*sin(100.000000) : -62.282974 , -97 123*sin(1000.000000) : 101.706184 , -25 123*sin(-0.010000) : -1.229980 , 0 123*sin(-0.100000) : -12.279510 , -9 123*sin(-1.000000) : -103.500931 , -100 123*sin(-10.000000) : 66.914597 , 59 123*sin(-100.000000) : 62.282974 , 99 123*sin(-1000.000000) : -101.706184 , 22 

why does this happen? is there a logical table error?

Answer 1

[change]

The code is experiencing problems because the angle increases by 360 degrees due to incorrect "modulo" arithmetic in the OP code. The product rad[i]*UCHAR_MAX/2/M_PI converted to an (8-bit) unsigned char , which is module 256, but the code is scaling tables and the code is UCHAR_MAX (255). The last point of this answer details the aspects of this, but it is clear that tables and code should use 256, not 255.

 unsigned char radToIndex; radToIndex=rad[i]*UCHAR_MAX/2/M_PI; // wrong scaling radToIndex=rad[i]*(UCHAR_MAX+1)/2/M_PI; // right 

Also, note that the OP code has undefined behavior when radToIndex == UCHAR_MAX , since this is an invalid index for int sini[UCHAR_MAX]; .

Using the above fix and 3 fixes below: table size 256, rounded index, rounding the sine value, use double for the results of creating the table in:

 123*sin(2.000000) : 111.843584 , 112 123*sin(4.000000) : -93.086707 , -93 123*sin(6.000000) : -34.368106 , -35 123*sin(8.000000) : 121.691064 , 121 123*sin(10.000000) : -66.914597 , -65 123*sin(-2.000000) : -111.843584 , -112 123*sin(-4.000000) : 93.086707 , 93 123*sin(-6.000000) : 34.368106 , 35 123*sin(-8.000000) : -121.691064 , -121 123*sin(-10.000000) : 66.914597 , 65 

---

The code also experiences double rounding or more valuable: double truncation.

radToIndex=rad[i]*UCHAR_MAX/2/M_PI; cropped in the direction 0. Thus, the index is made smaller, not the closest.

Creating a table sini[i]=sinf(i*2*M_PI/UCHAR_MAX)*MAGNIFICATION; also trimmed in the direction of 0. Thus, sini[] is made smaller, no closer than int .

To improve, simply round to the nearest using round() .

 sini[i] = (int) roundf(sinf(i*2*M_PI/UCHAR_MAX)*MAGNIFICATION); radToIndex= (int) round(rad[i]*UCHAR_MAX/2/M_PI); 

---

Typically, since float usually 24-bit precision and int likely sign of 31+, use double to create a table for additional improvements.

 sini[i] = (int) round(sin(i*2.0*M_PI/UCHAR_MAX)*MAGNIFICATION); 

---

In addition, we recommend using UCHAR_MAX + 1 See BAM :

Disabled at 1.

the index of the array becomes UCHAR_MAX, as it overflows, it automatically becomes 0

UCHAR_MAX does not overflow, UCHAR_MAX + 1 overflows and becomes 0. ( unsigned char math)

 int sini[UCHAR_MAX+1]; for (int i=0; i<(UCHAR_MAX+1); i++) { // Rather than `i*2*M_PI/UCHAR_MAX`, use sini[i]=sinf(i*2*M_PI/(UCHAR_MAX + 1))*MAGNIFICATION; 

Answer 2

Source of problem

It looks like you get errors from rounding floating point numbers and assign unsigned char floating point numbers.

The following program, adapted from your published code, demonstrates how the index begins to deviate even after rounding a floating point number.

 #include <limits.h> #include <math.h> int sini[UCHAR_MAX]; int cosi[UCHAR_MAX]; double angle[UCHAR_MAX]; #define MAGNIFICATION 256 #define SIN(i) sini[i]/MAGNIFICATION #define COS(i) cosi[i]/MAGNIFICATION void initTable() { double M_PI = 4.0*atan(1.0); for(int i=0;i<UCHAR_MAX;i++) { angle[i] = i*2*M_PI/UCHAR_MAX; sini[i]=sinf(angle[i])*MAGNIFICATION; cosi[i]=cosf(angle[i])*MAGNIFICATION; } } #include <stdio.h> void test3() { int radToIndexInt; unsigned char radToIndexChar; float radTemp; float rad[]={2.0f,4.0f,6.0f,8.0f,10.0f,-2.0f,-4.0f,-6.0f,-8.0f,-10.0f}; double M_PI = 4.0*atan(1.0); for(int i=0;i<sizeof(rad)/sizeof(float);i++) { radTemp = rad[i]*UCHAR_MAX/2/M_PI; radToIndexInt = round(radTemp); radToIndexInt %= UCHAR_MAX; if ( radToIndexInt < 0 ) { radToIndexInt += UCHAR_MAX; } radToIndexChar = round(radTemp); printf("radToIndexInt: %d, radToIndexChar: %d\n", radToIndexInt, radToIndexChar); } } int main() { initTable(); test3(); return 0; } 

The output of the specified program:

 radToIndexInt: 81, radToIndexChar: 81 radToIndexInt: 162, radToIndexChar: 162 radToIndexInt: 244, radToIndexChar: 244 radToIndexInt: 70, radToIndexChar: 69 radToIndexInt: 151, radToIndexChar: 150 radToIndexInt: 174, radToIndexChar: 175 radToIndexInt: 93, radToIndexChar: 94 radToIndexInt: 11, radToIndexChar: 12 radToIndexInt: 185, radToIndexChar: 187 radToIndexInt: 104, radToIndexChar: 106 

Decision

Using

  radToIndex=round(radTemp); radToIndex %= UCHAR_MAX; if ( radToIndex < 0 ) { radToIndex += UCHAR_MAX; } 

to calculate the index, I get very close answers:

Here the program, adapted again from your published code, demonstrates that using the above logic works.

 #include <limits.h> #include <math.h> int sini[UCHAR_MAX]; int cosi[UCHAR_MAX]; double angle[UCHAR_MAX]; #define MAGNIFICATION 256 #define SIN(i) sini[i]/MAGNIFICATION #define COS(i) cosi[i]/MAGNIFICATION void initTable() { double M_PI = 4.0*atan(1.0); for(int i=0;i<UCHAR_MAX;i++) { angle[i] = i*2*M_PI/UCHAR_MAX; sini[i]=sinf(angle[i])*MAGNIFICATION; cosi[i]=cosf(angle[i])*MAGNIFICATION; } } #include <stdio.h> void test2() { int radToIndex; float radTemp; int scalar=123; float rad[]={0.01f,0.1f,1.0f,10.0f,100.0f,1000.0f,-0.01f,-0.1f,-1.0f,-10.0f,-100.0f,-1000.0f}; double M_PI = 4.0*atan(1.0); printf("scalar=%d\n",scalar); for(int i=0;i<sizeof(rad)/sizeof(float);i++) { radTemp = rad[i]*UCHAR_MAX/2/M_PI; radToIndex=round(radTemp); radToIndex %= UCHAR_MAX; if ( radToIndex < 0 ) { radToIndex += UCHAR_MAX; } printf("%d*sin(%f) : %f , %d\n", scalar,rad[i],scalar*sinf(rad[i]),scalar*SIN(radToIndex)); } } void test1() { int radToIndex; float radTemp; int scalar=123; float rad[]={2.0f,4.0f,6.0f,8.0f,10.0f,-2.0f,-4.0f,-6.0f,-8.0f,-10.0f}; double M_PI = 4.0*atan(1.0); printf("scalar=%d\n",scalar); for(int i=0;i<sizeof(rad)/sizeof(float);i++) { radTemp = rad[i]*UCHAR_MAX/2/M_PI; radToIndex=round(radTemp); radToIndex %= UCHAR_MAX; if ( radToIndex < 0 ) { radToIndex += UCHAR_MAX; } printf("%d*sin(%f) : %f , %d\n", scalar,rad[i],scalar*sinf(rad[i]),scalar*SIN(radToIndex)); } } int main() { initTable(); test1(); test2(); return 0; } 

Output:

 scalar=123 123*sin(2.000000) : 111.843582 , 111 123*sin(4.000000) : -93.086708 , -92 123*sin(6.000000) : -34.368107 , -32 123*sin(8.000000) : 121.691063 , 121 123*sin(10.000000) : -66.914597 , -67 123*sin(-2.000000) : -111.843582 , -111 123*sin(-4.000000) : 93.086708 , 92 123*sin(-6.000000) : 34.368107 , 32 123*sin(-8.000000) : -121.691063 , -121 123*sin(-10.000000) : 66.914597 , 67 scalar=123 123*sin(0.010000) : 1.229980 , 0 123*sin(0.100000) : 12.279510 , 12 123*sin(1.000000) : 103.500931 , 103 123*sin(10.000000) : -66.914597 , -67 123*sin(100.000000) : -62.282974 , -63 123*sin(1000.000000) : 101.706184 , 102 123*sin(-0.010000) : -1.229980 , 0 123*sin(-0.100000) : -12.279510 , -12 123*sin(-1.000000) : -103.500931 , -103 123*sin(-10.000000) : 66.914597 , 67 123*sin(-100.000000) : 62.282974 , 63 123*sin(-1000.000000) : -101.706184 , -102