I am trying to implement the CORDIC algorithm for approximating a sine function in single precision, on architecture with no FPU. I compare the result obtained from my implementation with results obtained from standard C math functions. I tried implementing in two ways: 1) directly using floating-point operations, 2) converting the input to fixed-point and using integer based operations. I compare the results obtained from sinf(), sin() and sin() cast to float. The comparison is based on comparing the hexadecimal representations of the result with the expected from math functions.
In (1) the implementation uses double types, then the result is cast to float. My calculated values are always off by at least one hexadecimal digit, no matter how many iterations are done with CORDIC.
In (2), initially the input was mapped to a 32 bit integer. The error was the same as in (1). Only after, increasing the fixed point size to 64 bit (and the number of iterations to 64) the precision improved. But yet, there are ranges of input for which the algorithm is not precise. If I would increase the fixed point size to 128 bit (and the number of iterations to 128), it may be enough to obtain precise values but it is completely unpractical.
The algorithm in (1) is a modified version from the book 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;
}
The algorithm in (2) a is modified version found at 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;
}
The CORDIC table is similarly generated with bits=64.
The testing for (1) is done as below:
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);
}
The testing for (2) is done as below:
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);
}
Is there something that I don't take into account for CORDIC? Is it possible that CORDIC is completely not suitable and other methods should be considered?
I gave it a shot but as mentioned in the comments you can not expect exact bit match as math goniometrics is usually based on Chebyshev polynomials. Also you do not have defined the cordic_1K constant. After some search I managed to do this in 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);
}
//-------------------------------------------------------------------------
You can ignore the VCL stuff like AnsiString (or port it to your environment) its just for printing the results...
The code gave me this output:
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
As you can see the 64 bit tables produce much better match to math sin the error is up to 4 ulp (2^-24) for 32 bit and 2 ulp (2^-24) for 64 bit tables. As the mantissa of 32 bit float is 23+1 bits the result corresponds to 2 least significant bits of mantissa so its last hexadecimal digit is not matching ...
PS the atan table is:
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 };