Как се прилага Arctan?

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

Ето реализация на 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 обикновено са патентовани. За да изчислят arctan или други (обратни) тригонометрични функции, общите алгоритми следват процес от три стъпки:

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

Намаляването на аргумента обикновено се основава на добре известни тригонометрични идентичности, които могат да се търсят в различни стандартни препратки като MathWorld (http://mathworld.wolfram.com/InverseTangent.html). За изчисляването на arctan, често използваните идентичности са

• arctan (-x) = -arctan (x)
• 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ⁿ, което позволява използването на произволно тесен първичен интервал на приближение за сметка на допълнително съхранение на таблица. Това е класически програмен компромис между пространство и време.

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

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

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

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

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

По-долу е C код за изчисляване на arctan() с двойна точност, който демонстрира много от принципите и техниките на проектиране, споменати по-горе. На този бързо конструиран код липсва сложността на реализациите, посочени в други отговори, но трябва да осигури резултати с по-малко от 2 ulp грешки, което може да е достатъчно в различни контексти. Създадох персонализирано минимаксно приближение с проста реализация на алгоритъма на Remez, който използва 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) на ограничен брой интервали. Използвайте приближение на Чебишев за арктан от доста висок порядък за всеки интервал. (Това се прави офлайн и обикновено е източникът на всички магически числа, които виждате в тези реализации. Също така, човек може леко да затегне приближението на Чебишев, като използва алгоритъма за обмен на Remez, но не знам за случаи, в които това помага много .)
• Разберете в кой интервал е аргументът (използвайки ifs и други неща или просто трик с индексиране на таблица) и оценете серията на Чебишев на този интервал.

Няколко свойства са особено желани тук:

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

Не е напълно лесно да се оцени ред на Чебишев, така че тези две свойства да са валидни. Тук често се срещат трикове, при които две doubles се използват за представяне на различни части от една стойност. Тогава вероятно има някои казуси, които да покажат, че внедряването е монотонно. Също така, близо до нула, приближение на Тейлър до arctan вместо приближение на Чебишев --- търсите относителна грешка и оценяването на серията с помощта на правилото на Хорнър трябва да работи.

Ако търсите atan имплементация за четене, fdlibm изглежда по-малко гаден от този, който в момента е в glibc. Намаляването на аргумента изглежда се основава на самоличността на тригонометра tan(a+b) = (tan(a) + tan(b)) / (1 - tan(a) tan(b)), използвайки 0.5, 1 или 1.5 за tan(a) според случая.