Я пытаюсь реализовать алгоритм CORDIC для аппроксимации синусоидальной функции с одинарной точностью на архитектуре без FPU.Я сравниваю результат, полученный в моей реализации, с результатами, полученными от стандартных математических функций языка Си.Я попытался реализовать двумя способами: 1) напрямую используя операции с плавающей точкой, 2) преобразовать входные данные в операции с фиксированной точкой и использовать целочисленные операции.Я сравниваю результаты, полученные от sinf (), sin () и sin (), приведенных к float.Сравнение основано на сравнении шестнадцатеричных представлений результата с ожидаемыми математическими функциями.
В (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 аналогично генерируется сбиты = 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 совершенно не подходит, и следует рассмотреть другие методы?