c - 使用 CORDIC 的正弦计算结果不精确

Question

我正在尝试在没有 FPU 的架构上实现 CORDIC 算法以单精度逼近正弦函数。我将从我的实现中获得的结果与从标准 C 数学函数中获得的结果进行比较。我尝试以两种方式实现:1) 直接使用浮点运算，2) 将输入转换为定点并使用基于整数的运算。我比较了从 sinf()、sin() 和 sin() 转换为 float 获得的结果。比较基于将结果的十六进制表示与数学函数的预期值进行比较。

在 (1) 中，实现使用 double 类型，然后将结果转换为 float。无论使用 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)a中的算法是在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完全不适合，需要考虑其他方法？

Answer 1

我尝试了一下，但如评论中所述，您不能期望精确的位匹配，因为数学测角通常基于切比雪夫多项式。此外，您还没有定义 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);
    }
//-------------------------------------------------------------------------

您可以忽略 VCL 之类的东西，例如 AnsiString(或将其移植到您的环境中)，它只是为了打印结果...

代码给了我这个输出:

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) 和 2 ulp (2^-24) 用于 64 位表。由于 32 位 float 的尾数是 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 };