Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Why this sin cos look up table inaccurate when radian is large?

I want to create a sin cos look up table for optimization, by using an array index from 0 to UCHAR_MAX, so that 0 radian is index 0, pi/2 radian is 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 of using UCHAR_MAX as max is I want to make the good use of unsigned char overflow to simulates the radian thats varies from 0 to 2*pi only : for example, if the value of radian is 2*pi, the index of array becomes UCHAR_MAX, because it overflows, it automatically becomes 0 and no mod is required (if I use 0 to 360 as domain I may need to calculate index%360 every time). Then I test it with some radian values:

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

like the following:

#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),found the results starts go beyond the actual one when the magnitude of radian increases (when 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 test with another 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 increase when magnitude increases, so I am quite sure the table becomes inaccurate when radian is large. In sincos.h there is a value MAGNIFICATION to control the accuracy, I changed it from 256 to 4096, but it seems no much improvement:

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 would that happen? is there logical error of the table?

like image 595
ggrr Avatar asked Aug 31 '15 04:08

ggrr


People also ask

Should sin be in radians or degrees?

Most of the time we measure angles in degrees. For example, there are 360° in a full circle or one cycle of a sine wave, and sin(30°) = 0.5 and cos(90°) = 0. But it turns out that a more natural measure for angles, at least in mathematics, is in radians.

Can sin be in radians?

The value of sine 1 in radian is 0.8414709848.


1 Answers

[Edit]

Code experiences problems as the angle increases past 360 degrees due to wrong "modulo" arithmetic in OP's following code. The product rad[i]*UCHAR_MAX/2/M_PI is converted to an (8-bit) unsigned char which is a modulo 256, yet code is scaling tables and code by UCHAR_MAX (255). The last point of this answer details aspects of this, yet it is clear that tables and code should be 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

Further, note OP's code has undefined behavior when radToIndex == UCHAR_MAX as that is an invalid index to int sini[UCHAR_MAX];.

Using above fix and 3 below fixes: table size 256, round index, round the value of sine, use double for table creation results 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

Code is also experiencing double rounding or more preciously: double truncation.

radToIndex=rad[i]*UCHAR_MAX/2/M_PI; truncates toward 0. So the index is made smaller, not closest.

Table creation sini[i]=sinf(i*2*M_PI/UCHAR_MAX)*MAGNIFICATION; also truncates toward 0. So the sini[] is made smaller, not closest int.

To improve, simply round to nearest with round().

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

As a general note, since float is typically 24 bit precision and int likely 31+sign, use double for table creation for additional improvements.

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

Further, recommend using UCHAR_MAX + 1 See BAM:

Off by 1.

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

UCHAR_MAX is 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;
like image 182
chux - Reinstate Monica Avatar answered Sep 28 '22 10:09

chux - Reinstate Monica