logo标签列表

使用CORDIC计算正弦导致结果不准确

cmathnumericieee-754approximation
3
我正在尝试在没有FPU的体系结构上实现CORDIC算法以近似正弦函数的单精度计算。我将我的实现结果与标准的C数学函数结果进行了比较。我尝试了两种实现方式：1）直接使用浮点运算，2）将输入转换为定点并使用基于整数的运算。我比较了sinf()，sin()和转换为浮点数的sin()的结果。比较基于将结果的十六进制表示与来自数学函数的预期结果进行比较。
在(1)中，实现使用双精度类型，然后将结果转换为浮点数。无论使用CORDIC进行多少次迭代，我的计算值始终至少偏离一个十六进制数字。
在(2)中，最初将输入映射到32位整数。误差与(1)相同。只有在将固定点大小增加到64位（迭代次数增加到64）后，才提高了精度。但是，在某些输入范围内，该算法仍不准确。如果将固定点大小增加到128位（迭代次数增加到128），则可能足以获得准确值，但这完全不切实际。
(1)中的算法是从书籍https://www.jjj.de/fxt/fxtbook.pdf修改而来。
#include <math.h>
#include <stdio.h>
const double cordic_1K = 0.6072529350088812561694467525049282631123908521500897724;
double *cordic_ctab;

void make_cordic_ctab(ulong na)
{
    double s = 1.0;
    for (ulong k=0; k<na; ++k)
    {
        cordic_ctab[k] = atan(s);
        s *= 0.5;
    }
}


void cordic(int theta, double* s, double* c, int n)
{
    double x, y, z, v;
    double tx, ty, tz;
    double d;

    x = cordic_1K;
    y = 0;
    z = theta;
    v = 1.0;

    for (int k = 0; k < n; ++k) {
        d = (z >= 0 ? +1 : -1);
        tx = x - d * v * y;
        ty = y + d * v * x;
        tz = z - d * cordic_ctab[k];
        x = tx;
        y = ty;
        z = tz;
        v *= 0.5;
    }
    *c = x;
    *s = y;
}

(2)中的算法是修改版本，可在http://www.dcs.gla.ac.uk/~jhw/cordic/找到。

#include <math.h>
#include <stdio.h>
#define cordic_1K 0x26dd3b6a10d79600
#define CORDIC_NTAB 64

void cordic(long theta, long *s, long *c, int n) {
  long d, tx, ty, tz;
  long x = cordic_1K, y = 0, z = theta;
  n = (n > CORDIC_NTAB) ? CORDIC_NTAB : n;

  for (int k = 0; k < n; ++k) {
    d = z >= 0 ? 0 : -1;
    tx = x - (((y >> k) ^ d) - d);
    ty = y + (((x >> k) ^ d) - d);
    tz = z - ((cordic_ctab[k] ^ d) - d);
    x = tx;
    y = ty;
    z = tz;
  }

  *c = x;
  *s = y;
}

CORDIC表格同样是使用bits=64生成的。
对于（1）的测试如下进行：
int main(int argc, char **argv) {
  float angle;
  long s, c;
  int failed = 0;

  cordic_ctab = (double*)malloc(sizeof(double) * 64);
  make_cordic_ctab(64);

  for (int i = 0; i < step; i++) {
    angle = (i / step) * M_PI / 4;

    cordic(angle, &s, &c, 64);
    float result = s;
    float expected = sinf(angle);

    if (angle < pow(2, -27))
      result = angle;

    if (memcmp(&result, &expected, sizeof(float)) != 0) {
      failed += 1;
      printf("%e : %e\n", result, expected);
      printf("0x%x : 0x%x\n", *((unsigned int *)&result),
             *((unsigned int *)&expected));
      printf("\n");

    }
  }
  printf("failed:%d\n", failed);
}

对于(2)的测试如下：

int main(int argc, char **argv) {
  float angle;
  long s, c;
  int failed = 0;
  double mul = 4611686018427387904.000000;
  double step = 1000000000.0;

  for (int i = 0; i < step; i++) {
    angle = (i / step) * M_PI / 4;

    cordic((angle * mul), &s, &c, 64);
    float result = s / mul;
    float expected = sinf(angle);

    if (angle < pow(2, -27))
      result = angle;

    if (memcmp(&result, &expected, sizeof(float)) != 0) {
      failed += 1;
      printf("%e : %e\n", result, expected);
      printf("0x%x : 0x%x\n", *((unsigned int *)&result),
             *((unsigned int *)&expected));
      printf("\n");

    }
  }
  printf("failed:%d\n", failed);
}

我是否没有考虑到CORDIC的某些因素？是否有可能CORDIC完全不适用，需要考虑其他方法？

请展示您的代码、输入值以及算法输出结果，与标准库函数的输出结果进行比较。您所说的“单精度”是什么意思？您所说的“比较十六进制表示”是什么意思？如何获取这些信息？根据您处理器的能力，其他近似方法，例如泰勒级数，可能更快或更好。 - Bodo
@Bodo 我不同意你的观点。如果我假设另一个算法的结果是正确的，那么我的结果必须严格相同。否则，我就是在接近另一个近似算法的结果。如果我们讨论使用浮点数计算物理距离，那么是的，我需要比较差异并查看是否正确到某个足够的误差。 - Daniel
@MattTimmermans 为什么不呢？是CORDIC的分辨率不够好吗？还是我需要做更多的迭代？为什么使用泰勒级数可以使它成为位等效？ - Daniel
  1. CORDIC算法在整个范围内提供大致相同的精度，但浮点数和泰勒级数在接近零点时更加精确；
  2. 如果理想值非常接近两个浮点数的中间点，则即使是最小的误差也可能使其四舍五入到不同的方向。
 - Matt Timmermans
这正是我怀疑的，你证实了我的想法。谢谢。 - Daniel
显示剩余8条评论
1个回答
1
我试过了，但是正如评论中提到的，你不能期望精确的位匹配，因为数学三角函数通常基于Chebyshev多项式。此外，您还没有定义cordic_1K常量。经过一些搜索，我设法在C++/VCL中完成了这个任务。
//---------------------------------------------------------------------------
// IEEE 754 single masks
const DWORD _f32_sig    =0x80000000;    // sign
const DWORD _f32_exp    =0x7F800000;    // exponent
const DWORD _f32_exp_sig=0x40000000;    // exponent sign
const DWORD _f32_exp_bia=0x3F800000;    // exponent bias
const DWORD _f32_exp_lsb=0x00800000;    // exponent LSB
const DWORD _f32_exp_pos=        23;    // exponent LSB bit position
const DWORD _f32_man    =0x007FFFFF;    // mantisa
const DWORD _f32_man_msb=0x00400000;    // mantisa MSB
const DWORD _f32_man_bits=       23;    // mantisa bits
const float _f32_lsb     =  3.4e-38;    // abs min number
//---------------------------------------------------------------------------
float CORDIC32_atan[_f32_man_bits+1];
void f32_sincos(float &s,float &c,float a)
    {
    int k;
    float x,y=0.0,v=1.0,d,tx,ty,ta;
    x=0.6072529350088812561694; // cordic_1K
    for (k=0;k<=_f32_man_bits;k++)
        {
        d =(a>=0.0?+1.0:-1.0);
        tx=x-d*v*y;
        ty=y+d*v*x;
        ta=a-d*CORDIC32_atan[k];
        x=tx; y=ty; a=ta; v*=0.5;
        }
    c=x; s=y;
    }
//---------------------------------------------------------------------------
double CORDIC64_atan[_f32_man_bits+1];
void f64_sincos(float &s,float &c,double a)
    {
    int k;
    double x,y=0.0,v=1.0,d,tx,ty,ta;
    x=0.6072529350088812561694; // cordic_1K
    for (k=0;k<=_f32_man_bits;k++)
        {
        d =(a>=0.0?+1.0:-1.0);
        tx=x-d*v*y;
        ty=y+d*v*x;
        ta=a-d*CORDIC64_atan[k];
        x=tx; y=ty; a=ta; v*=0.5;
        }
    c=x; s=y;
    }
//---------------------------------------------------------------------------
//--- Builder: --------------------------------------------------------------
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    int i;
    float s0,c0,s1,c1,s2,c2,d32,d64,D32,D64,x;
    AnsiString txt="";

    // init CORDIC tables
    for (x=1.0,i=0;i<=_f32_man_bits;i++,x*=0.5)
        {
        CORDIC32_atan[i]=atan(x);
        CORDIC64_atan[i]=atan(double(x));
        }

    // 32 bit
    D32=0.0; D64=0.0;
    for (x=-0.5*M_PI;x<=+0.5*M_PI;x+=0.025)
        {
        s0=sin(x);
        c0=cos(x);
        f32_sincos(s1,c1,x); d32=fabs(s1-s0); if (D32<d32) D32=d32;
        f64_sincos(s2,c2,x); d64=fabs(s2-s0); if (D64<d64) D64=d64;
        if (d32+d64>1e-16)
            {
            txt+=AnsiString().sprintf("sin(%2.5f) == %2.5f != %2.5f != %2.5f  | %.10f %.10f\r\n",x,s0,s1,s2,d32,d64);
            f32_sincos(s0,c0,x);    // debug breakpoint
            f64_sincos(s2,c2,x);
            }
        }
    txt=AnsiString().sprintf("max err: %.10f %.10f\r\n",D32,D64)+txt;
    mm_log->Lines->Add(txt);
    }
//-------------------------------------------------------------------------

你可以忽略像 AnsiString 这样的VCL内容（或将其移植到您的环境中）仅用于打印结果... 代码给了我这个输出： 
max err: 0.0000002384 0.0000001192
sin(-1.54580) == -0.99969 != -0.99969 != -0.99969  | 0.0000000596 0.0000000000
sin(-1.52080) == -0.99875 != -0.99875 != -0.99875  | 0.0000001192 0.0000000000
sin(-1.49580) == -0.99719 != -0.99719 != -0.99719  | 0.0000001192 0.0000000000
sin(-1.44580) == -0.99220 != -0.99220 != -0.99220  | 0.0000000596 0.0000000000
sin(-1.42080) == -0.98877 != -0.98877 != -0.98877  | 0.0000000596 0.0000000596
sin(-1.39580) == -0.98473 != -0.98473 != -0.98473  | 0.0000001192 0.0000000596
sin(-1.37080) == -0.98007 != -0.98007 != -0.98007  | 0.0000000000 0.0000000596
sin(-1.34580) == -0.97479 != -0.97479 != -0.97479  | 0.0000000596 0.0000000000
sin(-1.27080) == -0.95534 != -0.95534 != -0.95534  | 0.0000001192 0.0000000000
sin(-1.24580) == -0.94765 != -0.94765 != -0.94765  | 0.0000000596 0.0000000596
sin(-1.22080) == -0.93937 != -0.93937 != -0.93937  | 0.0000000596 0.0000000000
sin(-1.19580) == -0.93051 != -0.93051 != -0.93051  | 0.0000000596 0.0000000596
sin(-1.17080) == -0.92106 != -0.92106 != -0.92106  | 0.0000000596 0.0000000000
sin(-1.14580) == -0.91104 != -0.91104 != -0.91104  | 0.0000001192 0.0000000596
sin(-1.12080) == -0.90045 != -0.90045 != -0.90045  | 0.0000001192 0.0000000000
sin(-1.07080) == -0.87758 != -0.87758 != -0.87758  | 0.0000001788 0.0000000596
sin(-1.04580) == -0.86532 != -0.86532 != -0.86532  | 0.0000001788 0.0000000596
sin(-1.02080) == -0.85252 != -0.85252 != -0.85252  | 0.0000001192 0.0000000596
sin(-0.99580) == -0.83919 != -0.83919 != -0.83919  | 0.0000000000 0.0000000596
sin(-0.97080) == -0.82534 != -0.82534 != -0.82534  | 0.0000001192 0.0000000596
sin(-0.94580) == -0.81096 != -0.81096 != -0.81096  | 0.0000000596 0.0000000596
sin(-0.92080) == -0.79608 != -0.79608 != -0.79608  | 0.0000000000 0.0000000596
sin(-0.89580) == -0.78071 != -0.78071 != -0.78071  | 0.0000001788 0.0000000596
sin(-0.87080) == -0.76484 != -0.76484 != -0.76484  | 0.0000000596 0.0000000000
sin(-0.84580) == -0.74850 != -0.74850 != -0.74850  | 0.0000000596 0.0000000000
sin(-0.82080) == -0.73169 != -0.73169 != -0.73169  | 0.0000001192 0.0000000000
sin(-0.79580) == -0.71442 != -0.71442 != -0.71442  | 0.0000000596 0.0000000000
sin(-0.77080) == -0.69671 != -0.69671 != -0.69671  | 0.0000000596 0.0000000596
sin(-0.74580) == -0.67856 != -0.67856 != -0.67856  | 0.0000000000 0.0000000596
sin(-0.72080) == -0.65998 != -0.65998 != -0.65998  | 0.0000001192 0.0000000596
sin(-0.69580) == -0.64100 != -0.64100 != -0.64100  | 0.0000000596 0.0000000000
sin(-0.67080) == -0.62161 != -0.62161 != -0.62161  | 0.0000001788 0.0000000596
sin(-0.64580) == -0.60184 != -0.60184 != -0.60184  | 0.0000000596 0.0000001192
sin(-0.62080) == -0.58168 != -0.58168 != -0.58168  | 0.0000000596 0.0000001192
sin(-0.59580) == -0.56117 != -0.56117 != -0.56117  | 0.0000000596 0.0000000596
sin(-0.57080) == -0.54030 != -0.54030 != -0.54030  | 0.0000001788 0.0000001192
sin(-0.54580) == -0.51910 != -0.51910 != -0.51910  | 0.0000001788 0.0000001192
sin(-0.52080) == -0.49757 != -0.49757 != -0.49757  | 0.0000000596 0.0000000894
sin(-0.49580) == -0.47573 != -0.47573 != -0.47573  | 0.0000000894 0.0000000596
sin(-0.47080) == -0.45360 != -0.45360 != -0.45360  | 0.0000000894 0.0000000298
sin(-0.44580) == -0.43118 != -0.43118 != -0.43118  | 0.0000000298 0.0000000298
sin(-0.42080) == -0.40849 != -0.40849 != -0.40849  | 0.0000000894 0.0000000596
sin(-0.39580) == -0.38554 != -0.38554 != -0.38554  | 0.0000001192 0.0000000596
sin(-0.37080) == -0.36236 != -0.36236 != -0.36236  | 0.0000000298 0.0000000000
sin(-0.34580) == -0.33895 != -0.33895 != -0.33895  | 0.0000000000 0.0000000298
sin(-0.32080) == -0.31532 != -0.31532 != -0.31532  | 0.0000000596 0.0000000000
sin(-0.29580) == -0.29150 != -0.29150 != -0.29150  | 0.0000000596 0.0000000596
sin(-0.27080) == -0.26750 != -0.26750 != -0.26750  | 0.0000000894 0.0000001192
sin(-0.24580) == -0.24333 != -0.24333 != -0.24333  | 0.0000000894 0.0000001192
sin(-0.22080) == -0.21901 != -0.21901 != -0.21901  | 0.0000000745 0.0000000894
sin(-0.19580) == -0.19455 != -0.19455 != -0.19455  | 0.0000000894 0.0000000596
sin(-0.17080) == -0.16997 != -0.16997 != -0.16997  | 0.0000001043 0.0000000894
sin(-0.14580) == -0.14528 != -0.14528 != -0.14528  | 0.0000000894 0.0000000894
sin(-0.12080) == -0.12050 != -0.12050 != -0.12050  | 0.0000000596 0.0000000671
sin(-0.09580) == -0.09565 != -0.09565 != -0.09565  | 0.0000000522 0.0000000522
sin(-0.07080) == -0.07074 != -0.07074 != -0.07074  | 0.0000000075 0.0000000224
sin(-0.04580) == -0.04578 != -0.04578 != -0.04578  | 0.0000000447 0.0000000335
sin(-0.02080) == -0.02080 != -0.02080 != -0.02080  | 0.0000000596 0.0000000577
sin(0.00420) == 0.00420 != 0.00420 != 0.00420  | 0.0000000545 0.0000000549
sin(0.02920) == 0.02920 != 0.02920 != 0.02920  | 0.0000000447 0.0000000410
sin(0.05420) == 0.05418 != 0.05418 != 0.05418  | 0.0000000149 0.0000000186
sin(0.07920) == 0.07912 != 0.07912 != 0.07912  | 0.0000000224 0.0000000373
sin(0.10420) == 0.10401 != 0.10401 != 0.10401  | 0.0000000820 0.0000000745
sin(0.12920) == 0.12884 != 0.12884 != 0.12884  | 0.0000001043 0.0000000894
sin(0.15420) == 0.15359 != 0.15359 != 0.15359  | 0.0000001043 0.0000001043
sin(0.17920) == 0.17825 != 0.17825 != 0.17825  | 0.0000000447 0.0000000745
sin(0.20420) == 0.20279 != 0.20279 != 0.20279  | 0.0000000596 0.0000000745
sin(0.22920) == 0.22720 != 0.22720 != 0.22720  | 0.0000001043 0.0000001043
sin(0.25420) == 0.25147 != 0.25147 != 0.25147  | 0.0000001192 0.0000000894
sin(0.27920) == 0.27559 != 0.27559 != 0.27559  | 0.0000000000 0.0000000596
sin(0.30420) == 0.29953 != 0.29953 != 0.29953  | 0.0000000596 0.0000000298
sin(0.32920) == 0.32329 != 0.32329 != 0.32329  | 0.0000000596 0.0000000596
sin(0.35420) == 0.34684 != 0.34684 != 0.34684  | 0.0000000298 0.0000000298
sin(0.37920) == 0.37018 != 0.37018 != 0.37018  | 0.0000000298 0.0000000298
sin(0.40420) == 0.39329 != 0.39329 != 0.39329  | 0.0000001788 0.0000000596
sin(0.42920) == 0.41615 != 0.41615 != 0.41615  | 0.0000000894 0.0000000894
sin(0.45420) == 0.43875 != 0.43875 != 0.43875  | 0.0000000298 0.0000000000
sin(0.47920) == 0.46107 != 0.46107 != 0.46107  | 0.0000000596 0.0000000298
sin(0.50420) == 0.48311 != 0.48311 != 0.48311  | 0.0000000596 0.0000000596
sin(0.52920) == 0.50485 != 0.50485 != 0.50485  | 0.0000001788 0.0000000596
sin(0.55420) == 0.52627 != 0.52627 != 0.52627  | 0.0000002384 0.0000001192
sin(0.57920) == 0.54736 != 0.54736 != 0.54736  | 0.0000001192 0.0000000596
sin(0.60420) == 0.56811 != 0.56811 != 0.56811  | 0.0000000596 0.0000000596
sin(0.62920) == 0.58850 != 0.58850 != 0.58850  | 0.0000000596 0.0000001192
sin(0.65420) == 0.60853 != 0.60853 != 0.60853  | 0.0000001192 0.0000001192
sin(0.67920) == 0.62817 != 0.62817 != 0.62817  | 0.0000000596 0.0000001192
sin(0.70420) == 0.64743 != 0.64743 != 0.64743  | 0.0000000596 0.0000000000
sin(0.72920) == 0.66628 != 0.66628 != 0.66628  | 0.0000000596 0.0000000000
sin(0.75420) == 0.68471 != 0.68471 != 0.68471  | 0.0000000596 0.0000000000
sin(0.77920) == 0.70271 != 0.70271 != 0.70271  | 0.0000000596 0.0000000000
sin(0.82920) == 0.73739 != 0.73739 != 0.73739  | 0.0000000596 0.0000000000
sin(0.85420) == 0.75405 != 0.75405 != 0.75405  | 0.0000001192 0.0000000000
sin(0.87920) == 0.77023 != 0.77023 != 0.77023  | 0.0000001192 0.0000000000
sin(0.90420) == 0.78593 != 0.78593 != 0.78593  | 0.0000000596 0.0000000596
sin(0.92920) == 0.80114 != 0.80114 != 0.80114  | 0.0000000596 0.0000001192
sin(0.95420) == 0.81585 != 0.81585 != 0.81585  | 0.0000001788 0.0000000596
sin(0.97920) == 0.83005 != 0.83005 != 0.83005  | 0.0000000000 0.0000000596
sin(1.00420) == 0.84373 != 0.84373 != 0.84373  | 0.0000001788 0.0000000000
sin(1.02920) == 0.85689 != 0.85689 != 0.85689  | 0.0000000596 0.0000000000
sin(1.05420) == 0.86951 != 0.86951 != 0.86951  | 0.0000001192 0.0000000000
sin(1.12920) == 0.90407 != 0.90407 != 0.90407  | 0.0000000596 0.0000000000
sin(1.15420) == 0.91447 != 0.91447 != 0.91447  | 0.0000000596 0.0000000596
sin(1.17920) == 0.92430 != 0.92430 != 0.92430  | 0.0000001788 0.0000000596
sin(1.20420) == 0.93355 != 0.93355 != 0.93355  | 0.0000000596 0.0000000000
sin(1.25420) == 0.95030 != 0.95030 != 0.95030  | 0.0000000596 0.0000000596
sin(1.27920) == 0.95779 != 0.95779 != 0.95779  | 0.0000001192 0.0000000596
sin(1.30420) == 0.96467 != 0.96467 != 0.96467  | 0.0000000596 0.0000000000
sin(1.35420) == 0.97664 != 0.97664 != 0.97663  | 0.0000000596 0.0000000596
sin(1.45420) == 0.99321 != 0.99321 != 0.99321  | 0.0000000596 0.0000000000
sin(1.47920) == 0.99581 != 0.99581 != 0.99581  | 0.0000001192 0.0000000000
sin(1.50420) == 0.99778 != 0.99778 != 0.99778  | 0.0000000596 0.0000000000
sin(1.52920) == 0.99914 != 0.99914 != 0.99914  | 0.0000000596 0.0000000000
sin(1.55420) == 0.99986 != 0.99986 != 0.99986  | 0.0000000596 0.0000000000

正如您所看到的，64位表对数学中的sin提供了更好的匹配，误差为32位的4 ulp（2^-24）和64位的2 ulp（2^-24）。由于32位浮点数的尾数是23+1位，因此结果对应于尾数的2个最低有效位，因此其最后一个十六进制数字不匹配...

PS atan表格如下：

CORDIC64_atan[24]={ 0.785398163397448, 0.463647609000806, 0.244978663126864, 0.124354994546761, 0.0624188099959574, 0.0312398334302683, 0.0156237286204768, 0.00781234106010111, 0.00390623013196697, 0.00195312251647882, 0.00097656218955932, 0.000488281211194898, 0.000244140620149362, 0.00012207031189367, 6.10351561742088E-05, 3.05175781155261E-05, 1.52587890613158E-05, 7.62939453110197E-06, 3.8146972656065E-06, 1.90734863281019E-06, 9.53674316405961E-07, 4.76837158203089E-07, 2.38418579101558E-07, 1.19209289550781E-07 };
谢谢！这意味着需要更多的迭代才能匹配数学中的sine函数吗？ - Daniel
@Daniel 不是的，每次迭代都代表尾数的一位，因此如果你进行的迭代次数超过了尾数位数，结果不会受到舍入的影响。差异主要是由于计算方式不同造成的。你知道数学中的正弦余弦也是有舍入误差的... - Spektre
对，我的意思是通过提高计算精度来增加迭代次数。只是为了完全理解：如果我可以在X位上计算CORDIC，其中X是相关范围内可能指数的数量乘以23位有效数字，那么我就可以得到浮点数的正确结果？ - Daniel
如果我没记错，数学中的 sin 和 cos 函数有 +/-1 ulp 的误差……因此，与 CORDIC 相比，最多可能存在 2 ulp 的误差差异，这可能只是两种变体都进行了相反的舍入，每个变体的舍入误差为 1 ulp……因此，64 位变体可能具有与您正在比较的 sin 和 cos 函数相同的精度…… - Spektre
1
@Spektre：关于“我记得数学中的sin和cos有+/-1 ulp误差”的问题：数学库函数的精度取决于实现。尽管我认为CRlibm覆盖了正弦和余弦，但没有人使用已知有界运行时间实现了所有的数学库函数。商业实现可能是3-4 ULP，例如。使用双精度sin作为单精度sinf的参考，可能会给出一个非常接近½ ULP值（使用单精度ULP）。 - Eric Postpischil
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接