使用 CORDIC 的正弦计算结果不精确
sine computation with CORDIC results in not precise values
我正在尝试在没有 FPU 的架构上实现 CORDIC 算法以单精度逼近正弦函数。我将从我的实现中获得的结果与从标准 C 数学函数中获得的结果进行比较。我尝试以两种方式实现：1) 直接使用浮点运算，2) 将输入转换为定点并使用基于整数的运算。
我比较了从 sinf()、sin() 和 sin() 转换为浮点数获得的结果。比较基于将结果的十六进制表示与数学函数的预期值进行比较。
(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 table 类似地生成，位=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完全不suitable，应该考虑其他方法？
我试了一下，但正如评论中提到的，您不能期望精确的位匹配，因为数学测角通常基于切比雪夫多项式。此外，您还没有定义 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 位 table 与数学 sin 的匹配更好，32 位的误差高达 4 ulp (2^-24) 并且 2 ulp (2^-24) 对于 64 位 tables。由于 32 位浮点数的尾数是 23+1 位，结果对应于尾数的 2 个最低有效位，因此它的最后一个十六进制数字不匹配 ...
PS atan table 是：
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 };
我正在尝试在没有 FPU 的架构上实现 CORDIC 算法以单精度逼近正弦函数。我将从我的实现中获得的结果与从标准 C 数学函数中获得的结果进行比较。我尝试以两种方式实现：1) 直接使用浮点运算，2) 将输入转换为定点并使用基于整数的运算。 我比较了从 sinf()、sin() 和 sin() 转换为浮点数获得的结果。比较基于将结果的十六进制表示与数学函数的预期值进行比较。
(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 table 类似地生成，位=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完全不suitable，应该考虑其他方法？
我试了一下，但正如评论中提到的，您不能期望精确的位匹配，因为数学测角通常基于切比雪夫多项式。此外，您还没有定义 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 位 table 与数学 sin 的匹配更好，32 位的误差高达 4 ulp (2^-24) 并且 2 ulp (2^-24) 对于 64 位 tables。由于 32 位浮点数的尾数是 23+1 位，结果对应于尾数的 2 个最低有效位，因此它的最后一个十六进制数字不匹配 ...
PS atan table 是：
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 };