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?
