Реализация Арксинус CORDIC не

Арксин "одиночного поворота" сильно ошибается, когда аргумент просто больше начального значения 'x', где это магический масштабирующий коэффициент - 1/An ~= 0.607252935 ~= 0x26DD3B6A.

Это происходит потому, что для всех аргументов > 0 Первый Шаг Всегда имеет y = 0

• Если arg

• Если arg > 1 / An, то d = +1, и этот шаг отодвигается все дальше от правильного ответа, а для диапазона значений чуть больше 1 / An все последующие шаги имеют d = -1, но не могут исправить результат :-(

Я нашел:

 arg = 0.607 (ie 0x26D91687), relative error 7.139E-09 -- OK    
 arg = 0.608 (ie 0x26E978D5), relative error 1.550E-01 -- APALLING !!
 arg = 0.685 (ie 0x2BD70A3D), relative error 2.667E-04 -- BAD !!
 arg = 0.686 (ie 0x2BE76C8B), relative error 1.232E-09 -- OK, again

Описания метода предупреждают об abs (arg) >= 0.98 (или около того), и я обнаружил, что где-то после 0.986 процесс не сходится и относительная ошибка скачет до ~5E-02 и достигает 1E-01 (!!) при arg=1 :- (

Как и вы, я также обнаружил, что для 0,303

Итак... единственный поворот CORDIC для arcsine выглядит для меня ерундой :- (

---

Добавлено позже... когда я еще внимательнее присмотрелся к единичному вращающемуся КОРДИКУ, то обнаружил еще много небольших областей, где относительная погрешность велика...

... поэтому я бы не касайтесь этого как метода вообще... это не просто мусор, этобесполезно .

---

Кстати: я очень рекомендую "руководство по программному обеспечению для элементарных функций", William Cody and William Waite, Prentice-Hall, 1980. Методы вычисления функций уже не так интересны (но требуется тщательное практическое обсуждение соответствующих сокращений диапазона). Однако для каждой функции они дают хорошую процедуру проверки.

Дополнительный источник, который я связал в конце вопроса, по-видимому, содержит решение.

Предлагаемый код можно свести к следующему:

#define M_PI_2_32    1.57079632F
#define SQRT2_2      7.071067811865476e-001F /* sin(45°) = cos(45°) = sqrt(2)/2 */

FLOAT32 angles[] = {
    7.8539816339744830962E-01F, 4.6364760900080611621E-01F, 2.4497866312686415417E-01F, 1.2435499454676143503E-01F,
    6.2418809995957348474E-02F, 3.1239833430268276254E-02F, 1.5623728620476830803E-02F, 7.8123410601011112965E-03F,
    3.9062301319669718276E-03F, 1.9531225164788186851E-03F, 9.7656218955931943040E-04F, 4.8828121119489827547E-04F,
    2.4414062014936176402E-04F, 1.2207031189367020424E-04F, 6.1035156174208775022E-05F, 3.0517578115526096862E-05F,
    1.5258789061315762107E-05F, 7.6293945311019702634E-06F, 3.8146972656064962829E-06F, 1.9073486328101870354E-06F,
    9.5367431640596087942E-07F, 4.7683715820308885993E-07F, 2.3841857910155798249E-07F, 1.1920928955078068531E-07F,
    5.9604644775390554414E-08F, 2.9802322387695303677E-08F, 1.4901161193847655147E-08F, 7.4505805969238279871E-09F,
    3.7252902984619140453E-09F, 1.8626451492309570291E-09F, 9.3132257461547851536E-10F, 4.6566128730773925778E-10F};

FLOAT32 arcsin_cordic(FLOAT32 t)
{            
    INT32 i;
    INT32 j;
    INT32 flip;
    FLOAT32 poweroftwo;
    FLOAT32 sigma;
    FLOAT32 sign_or;
    FLOAT32 theta;
    FLOAT32 x1;
    FLOAT32 x2;
    FLOAT32 y1;
    FLOAT32 y2;

    flip       = 0; 
    theta      = 0.0F;
    x1         = 1.0F;
    y1         = 0.0F;
    poweroftwo = 1.0F;

    /* If the angle is small, use the small angle approximation */
    if ((t >= -0.002F) && (t <= 0.002F))
    {
        return t;
    }

    if (t >= 0.0F) 
    {
        sign_or = 1.0F;
    }
    else
    {
        sign_or = -1.0F;
    }

    /* The inv_sqrt() is the famous Fast Inverse Square Root from the Quake 3 engine
       here used with 3 (!!) Newton iterations */
    if ((t >= SQRT2_2) || (t <= -SQRT2_2))
    {
        t =  1.0F/inv_sqrt(1-t*t);
        flip = 1;
    }

    if (t>=0.0F) 
    {
        sign_or = 1.0F;
    }
    else
    {
        sign_or = -1.0F;
    }

    for ( j = 0; j < 32; j++ ) 
    {
        if (y1 > t)
        {
            sigma = -1.0F;
        }
        else
        {
            sigma = 1.0F;
        }

        /* Here a double iteration is done */
        x2 =                       x1  - (sigma * poweroftwo * y1);
        y2 = (sigma * poweroftwo * x1) +                       y1;

        x1 =                       x2  - (sigma * poweroftwo * y2);
        y1 = (sigma * poweroftwo * x2) +                       y2;

        theta  += 2.0F * sigma * angles[j];

        t *= (1.0F + poweroftwo * poweroftwo);

        poweroftwo *= 0.5F;
    }

    /* Remove bias */
    theta -= sign_or*4.85E-8F;

    if (flip)
    {
        theta = sign_or*(M_PI_2_32-theta);
    }

    return theta;
}

Следует отметить следующее:

• это КОРДИЧЕСКАЯ реализация "двойной итерации".
Таким образом, таблица
• отличается по конструкции от старой таблицы.
• и вычисление выполняется в нотации с плавающей запятой, это приведет к значительному увеличению времени вычисления на целевое оборудование.
• на выходе присутствует небольшое смещение, удаляемое через проход theta -= sign_or*4.85E-8F;.

На следующем рисунке показаны абсолютные (слева) и относительные ошибки (справа) старой реализации (вверху) по сравнению с реализацией, содержащейся в этом ответе (внизу). Относительная погрешность получается только путем деления выхода CORDIC на выход встроенного математического аппарата.реализация сек. Она строится вокруг 1, а не 0 по этой причине.

Пик относительная погрешность (при не делении на ноль) равна 1.0728836e-006.

Средняя относительная погрешность равна 2.0253509e-007 (почти в соответствии с 32-битной точностью).