logo标签列表

为什么当弧度很大时，这个正弦余弦查找表不准确？

carraysmathtrigonometry
5

我希望创建一个正弦余弦查找表以进行优化，通过使用从0到UCHAR_MAX的数组索引，使得0弧度是索引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;
    }
}

使用UCHAR_MAX作为最大值的原因是我想充分利用无符号字符溢出来模拟弧度，它只变化从0到2*pi ：例如，如果弧度的值是2*pi，数组的索引就变成了UCHAR_MAX，因为它溢出了，它自动变成了0，不需要进行模运算（如果我使用0到360作为域，每次都需要计算index%360）。然后我使用一些弧度值进行测试：
float rad[]={2.0f,4.0f,6.0f,8.0f,10.0f,-2.0f,-4.0f,-6.0f,-8.0f,-10.0f};

像下面这样：

#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;
}

我使用 123*sin(radian) 测试了这个表格，发现当弧度的大小增加（当弧度为10或-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

并使用其他数据进行测试：

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};

输出：

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

当幅值增加时，误差也会增加，所以我非常确定当弧度很大时，表格会变得不准确。 在 sincos.h 中有一个名为 MAGNIFICATION 的值，用于控制精度。 我将其从 256 更改为 4096，但似乎没有多少改善：

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

为什么会发生这种情况？表格中是否存在逻辑错误？
你如何衡量准确性？ - Iman Nia
2个回答
5

[编辑]

由于OP以下代码中的错误“模运算”，当角度超过360度时，代码会出现问题。 乘积rad[i]*UCHAR_MAX/2/M_PI被转换为（8位）unsigned char，它是模256的，但是代码正在通过UCHAR_MAX（255）缩放表格和代码。 本答案的最后一点详细说明了这一点，但很明显，应该使用256而不是255来使用表格和代码。

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

此外，请注意当radToIndex == UCHAR_MAX时，OP的代码具有未定义的行为，因为那是一个无效的索引到int sini[UCHAR_MAX];。
使用上述修复和下面三个修复：表大小256，将索引舍入，将正弦值四舍五入，使用双精度浮点数创建表格，结果如下：
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

代码也会遇到“双重舍入”或者更准确地说是“双重截断”，double rounding。

radToIndex=rad[i]*UCHAR_MAX/2/M_PI; 截断至0，所以索引变小了，而不是最接近的。

表格创建 sini[i]=sinf(i*2*M_PI/UCHAR_MAX)*MAGNIFICATION; 也会截断至0。所以，sini[] 变小了，而不是最接近的 int。

要改进，只需使用 round() 进行四舍五入即可。

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

作为一般性的提示，由于float通常只有24位精度，而int则有31位加符号位，因此在表格创建时使用double可以获得额外的改进。
sini[i] = (int) round(sin(i*2.0*M_PI/UCHAR_MAX)*MAGNIFICATION);

另外，建议使用UCHAR_MAX + 1。参见BAM：

有一个偏差。

数组的索引变成了UCHAR_MAX，因为它溢出了，所以它自动变成了0。

UCHAR_MAX没有溢出，UCHAR_MAX + 1 溢出并变成了0。（unsigned char数学）

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;
在创建表格时使用double和sin()函数，当MAGNIFICATION > 1e23时非常有用。 - chux - Reinstate Monica
0

问题来源

看起来您正在遇到浮点数舍入和将浮点数赋值给unsigned char导致的错误。

以下程序是根据您发布的代码进行调整的，演示了即使在浮点数舍入后，索引也开始偏离的情况。

#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;
}

以上程序的输出结果为：
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

解决方案

通过使用

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

计算索引时，我得到了非常接近的答案：

这是一个程序，再次改编自您发布的代码，演示了使用上述逻辑的有效性。

#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;
}

输出：

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
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接