Как реализован арктан?

Тригонометрические функции имеют довольно уродливые реализации, которые являются хакерскими и делают много битов. Я думаю, что будет довольно сложно найти здесь кого-то, кто сможет объяснить алгоритм, который на самом деле используется.

Вот реализация atan2: https://sourceware.org/git/?p=glibc.git;a=blob;f=sysdeps/ieee754/dbl-64/e_atan2.c;h=a287ca6656b210c77367eec3c46d72f18476d61d;hb=HEAD

Изменить: на самом деле я нашел это: http://www.netlib.org/fdlibm/e_atan2.c за которым намного проще следить, но, вероятно, из-за этого он медленнее (?).

FPU делает все это в некоторых схемах, поэтому процессору не нужно выполнять всю эту работу.

Реализации инструкций FPATAN в процессорах x86 обычно являются проприетарными. Чтобы вычислить арктангенс или другие (обратные) тригонометрические функции, общие алгоритмы следуют трехэтапному процессу:

1. уменьшение аргумента для сопоставления полной входной области с узким интервалом
2. расчет базовой аппроксимации на узком интервале (интервал первичной аппроксимации)
3. расширение промежуточного результата на основе уменьшения аргумента для получения конечного результата

Сокращение аргументов обычно основано на хорошо известных тригонометрических тождествах, которые можно найти в различных стандартных справочниках, таких как MathWorld (http://mathworld.wolfram.com/InverseTangent.html). Для вычисления арктангенса обычно используются следующие тождества:

• арктангенс (-х) = -арктангенс (х)
• arctan (1/x) = 0,5 * pi - arctan(x) [x > 0]
• arctan (x) = arctan(c) + arctan((x - c) / (1 + x*c))

Обратите внимание, что последнее тождество поддается построению таблицы значений arctan(i/2ⁿ), i = 1...2ⁿ, которая позволяет использовать сколь угодно узкого интервала первичной аппроксимации за счет дополнительного табличного хранения. Это классический компромисс программирования между пространством и временем.

Аппроксимация на основном интервале обычно представляет собой минимаксную полиномиальную аппроксимацию достаточной степени. Рациональные приближения обычно неконкурентоспособны на современном оборудовании из-за высокой стоимости деления с плавающей запятой, а также страдают от дополнительной числовой ошибки из-за вычисления двух полиномов плюс ошибка, вносимая делением.

Коэффициенты для минимаксных полиномиальных аппроксимаций обычно вычисляются с использованием алгоритма Ремеза (http://en.wikipedia.org/wiki/Remez_algorithm). Такие инструменты, как Maple и Mathematica, имеют встроенные средства для вычисления таких приближений. Точность полиномиальных аппроксимаций можно повысить, убедившись, что все коэффициенты точно представляют машинные числа. Единственный известный мне инструмент со встроенной возможностью для этого — Sollya (http://sollya.gforge.inria.fr/), который предлагает функцию fpminimax().

При вычислении полиномов обычно используется схема Хорнера (http://en.wikipedia.org/wiki/Horner%27s_method), которая является эффективной и точной, или смесь схемы Эстрина (http://en.wikipedia.org/wiki/Estrin%27s_scheme) и Хорнера. Схема Эстрина позволяет превосходно использовать параллелизм на уровне команд, обеспечиваемый суперскалярными процессорами, с незначительным влиянием на общее количество инструкций и часто (но не всегда) благотворно влияет на точность.

Использование FMA (сложение с плавным умножением) повышает точность и производительность любой из схем оценки за счет уменьшения количества шагов округления и обеспечения некоторой защиты от вычитания. FMA встречается на многих процессорах, включая графические процессоры и последние процессоры x86. В стандартном C и стандартном C++ операция FMA представлена ​​как стандартная библиотечная функция fma(), однако ее необходимо эмулировать на платформах, которые не предлагают аппаратную поддержку, что замедляет ее работу на этих платформах.

С точки зрения программирования хотелось бы избежать риска ошибок преобразования при переводе констант с плавающей запятой, необходимых для аппроксимации и сокращения аргументов из текстового представления в машинное. Процедуры преобразования ASCII в числа с плавающей запятой печально известны наличием коварных ошибок (например, http://www.exploringbinary.com/php-hangs-on-numeric-value-2-2250738585072011e-308/). Один механизм, предлагаемый стандартным C (не C++, насколько я знаю, где он доступен только как проприетарное расширение), заключается в указании констант с плавающей запятой в виде шестнадцатеричных литералов, которые непосредственно выражают базовый битовый шаблон, эффективно избегая сложных преобразований.

Ниже приведен код C для вычисления arctan() с двойной точностью, который демонстрирует многие из принципов и методов проектирования, упомянутых выше. Этому быстро сконструированному коду не хватает сложности реализаций, указанных в других ответах, но он должен давать результаты с ошибкой менее 2 ulps, чего может быть достаточно в различных контекстах. Я создал собственное минимаксное приближение с простой реализацией алгоритма Ремеза, который использовал 1024-битную арифметику с плавающей запятой для всех промежуточных шагов. Я ожидаю, что использование Sollya или аналогичных инструментов приведет к численно превосходным приближениям.

double my_atan (double x)
{
    double a, z, p, r, s, q, o;
    /* argument reduction: 
       arctan (-x) = -arctan(x); 
       arctan (1/x) = 1/2 * pi - arctan (x), when x > 0
    */
    z = fabs (x);
    a = (z > 1.0) ? 1.0 / z : z;
    /* evaluate minimax polynomial approximation */
    s = a * a; // a**2
    q = s * s; // a**4
    o = q * q; // a**8
    /* use Estrin's scheme for low-order terms */
    p = fma (fma (fma (-0x1.53e1d2a25ff34p-16, s, 0x1.d3b63dbb65af4p-13), q,
                  fma (-0x1.312788dde0801p-10, s, 0x1.f9690c82492dbp-9)), o,
             fma (fma (-0x1.2cf5aabc7cef3p-7, s, 0x1.162b0b2a3bfcep-6), q, 
                  fma (-0x1.a7256feb6fc5cp-6, s, 0x1.171560ce4a483p-5)));
    /* use Horner's scheme for high-order terms */
    p = fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (fma (p, s,
        -0x1.4f44d841450e1p-5), s,
         0x1.7ee3d3f36bb94p-5), s, 
        -0x1.ad32ae04a9fd1p-5), s,  
         0x1.e17813d66954fp-5), s, 
        -0x1.11089ca9a5bcdp-4), s,  
         0x1.3b12b2db51738p-4), s,
        -0x1.745d022f8dc5cp-4), s,
         0x1.c71c709dfe927p-4), s,
        -0x1.2492491fa1744p-3), s,
         0x1.99999999840d2p-3), s,
        -0x1.555555555544cp-2) * s, a, a);
    /* back substitution based on argument reduction */
    r = (z > 1.0) ? (0x1.921fb54442d18p+0 - p) : p;
    return copysign (r, x);
}

Резюме: Тяжело. Кроме того, Эрик Постпишил и Стивен Кэнон, которые иногда околачиваются вокруг SO, очень хороши в этом.

Обычный подход для многих специальных функций выглядит следующим образом:

• Обрабатывайте NaN, бесконечности и нули со знаком как особые случаи.
• Если число настолько велико, что результат округляется до M_PI, вернуть M_PI. Назовите этот порог M.
• Если есть какая-либо идентичность сокращения аргумента, используйте ее, чтобы привести аргумент в более удобный диапазон. (Это может быть сложно: для sin и cos это означает, что вы выбираете кратное точному значение 2pi, чтобы попасть в правильный диапазон.)
• Разбейте [0,M) на конечное число интервалов. Используйте аппроксимацию Чебышева для арктангенса довольно высокого порядка на каждом интервале. (Это делается в автономном режиме, и обычно это источник всех магических чисел, которые вы видите в этих реализациях. Кроме того, можно немного подтянуть аппроксимацию Чебышева, используя алгоритм обмена Ремеза, но я не знаю ни одного случая, когда это сильно помогает. .)
• Выясните, в каком интервале находится аргумент (используя ifs и прочее, или просто трюк с индексацией таблицы), и оцените ряд Чебышева на этом интервале.

Здесь особенно желательны несколько свойств:

• Реализация arctan должна быть монотонной; то есть если x < y, то arctan(x) <= arctan(y).
• Реализация arctan всегда должна возвращать ответ в пределах 1 ulp от правильного ответа. Обратите внимание, что это относительная ошибка.

Не совсем просто вычислить ряд Чебышева так, чтобы выполнялись эти два свойства. Здесь распространены приемы, когда два double используются для представления разных частей одного значения. Тогда, вероятно, есть некоторые примеры, чтобы показать, что реализация монотонна. Кроме того, около нуля приближение Тейлора к arctan вместо приближения Чебышева --- вы после границы относительной ошибки и оценка ряда с использованием правила Хорнера должна работать.

Если вы ищете реализацию atan для чтения, fdlibm кажется менее неприятной, чем та, что в настоящее время находится в glibc. Сокращение аргумента, по-видимому, основано на идентификаторе триггера tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a) tan(b)) с использованием 0.5, 1 или 1.5 вместо tan(a) в зависимости от ситуации.