Недавно я реализовал библиотеку функций CORDIC, чтобы уменьшить требуемую вычислительную мощность (мой проект основан на PowerPC и очень строгий по своим спецификациям времени выполнения). Язык - это ANSI-C.
Другие функции (sin/cos/atan) работают с точностью до пределов точности как в 32, так и в 64-битных реализациях.К сожалению, функция asin() систематически не работает для определенных входов.
В целях тестирования я применил файл .h, который будет использоваться в симуляционной S-функции. (Это только для моего удобства, вы можете скомпилировать следующее как самостоятельное .exe с минимальными изменениями)
Примечание. Я заставил 32 итерации, потому что я работаю с 32-битной точностью и требуется максимально возможная точность.
Cordic.h:
#include <stdio.h>
#include <stdlib.h>
#define FLOAT32 float
#define INT32 signed long int
#define BIT_XOR ^
#define CORDIC_1K_32 0x26DD3B6A
#define MUL_32 1073741824.0F /*needed to scale float -> int*/
#define INV_MUL_32 9.313225746E-10F /*needed to scale int -> float*/
INT32 CORDIC_CTAB_32 [] = {0x3243f6a8, 0x1dac6705, 0x0fadbafc, 0x07f56ea6, 0x03feab76, 0x01ffd55b, 0x00fffaaa, 0x007fff55,
0x003fffea, 0x001ffffd, 0x000fffff, 0x0007ffff, 0x0003ffff, 0x0001ffff, 0x0000ffff, 0x00007fff,
0x00003fff, 0x00001fff, 0x00000fff, 0x000007ff, 0x000003ff, 0x000001ff, 0x000000ff, 0x0000007f,
0x0000003f, 0x0000001f, 0x0000000f, 0x00000008, 0x00000004, 0x00000002, 0x00000001, 0x00000000};
/* CORDIC Arcsine Core: vectoring mode */
INT32 CORDIC_asin(INT32 arc_in)
{
INT32 k;
INT32 d;
INT32 tx;
INT32 ty;
INT32 x;
INT32 y;
INT32 z;
x=CORDIC_1K_32;
y=0;
z=0;
for (k=0; k<32; ++k)
{
d = (arc_in - y)>>(31);
tx = x - (((y>>k) BIT_XOR d) - d);
ty = y + (((x>>k) BIT_XOR d) - d);
z += ((CORDIC_CTAB_32[k] BIT_XOR d) - d);
x = tx;
y = ty;
}
return z;
}
/* Wrapper function for scaling in-out of cordic core*/
FLOAT32 asin_wrap(FLOAT32 arc)
{
return ((FLOAT32)(CORDIC_asin((INT32)(arc*MUL_32))*INV_MUL_32));
}
Это можно вызвать аналогично:
#include "Cordic.h"
#include "math.h"
void main()
{
y1 = asin_wrap(value_32); /*my implementation*/
y2 = asinf(value_32); /*standard math.h for comparison*/
}
Результаты показаны так:
В верхнем левом углу находится вкладка [-1;1] более чем на 2000 шагов (0.001), в нижней части слева выводится моя функция, внизу справа - стандартный вывод, а верхний правый - разность двух выходов.
Немедленно увидеть, что ошибка не соответствует 32-битной точности.
Я проанализировал выполненные шаги (и промежуточные результаты) по моему коду, и мне кажется, что в какой-то момент значение y "достаточно близко" к начальному значению arc_in, и что может быть связанными с бит-сдвигом, приводит к расходимости решения.
Мои вопросы:
- Я в недоумении, это ошибка, присущая реализации CORDIC, или я допустил ошибку в реализации? Я ожидал снижения точности вблизи экстремумов, но эти всплески в середине довольно неожиданны. (наиболее заметные из них находятся за пределами
+/- 0.6, но даже удалены они больше при меньших значениях, хотя и не столь выражены) - Если это какая-то часть реализации CORDIC, существуют ли известные способы обхода?
EDIT:
-
Так как некоторые комментарии упоминают об этом, да, я тестировал определение
INT32, даже записывая#define INT32 int32_Tне изменяет результаты по малейшей сумме. -
Время вычисления на целевом оборудовании было измерено сотнями повторений блока из 10.000 итераций функции со случайным входом в диапазоне допустимости. Наблюдаемые средние результаты (для одного вызова функции) заключаются в следующем:
math.h asinf() 100.00 microseconds CORDIC asin() 5.15 microseconds
(по-видимому, предыдущий тест был неисправен, новый кросс-тест получил не лучше, чем в среднем 100 микросекунд в диапазоне допустимости) -
Я, по-видимому, нашел лучшую реализацию. Его можно загрузить в версии matlab здесь и в C здесь. Я буду более подробно анализировать его внутренние работы и сообщать позже.
