Die approximation der inversen trigonometrischen Funktionen

In einem fixed-point-Umgebung (S15.16) ich habe erfolgreich verwendet den CORDIC-Algorithmus (siehe Wikipedia für eine Allgemeine Beschreibung) zu berechnen atan2(y,x), dann aus der asin() und acos() aus, die mit bekannten funktionalen Identitäten, die mit der Quadratwurzel:

asin(x) = atan2 (x, sqrt ((1.0 + x) * (1.0 - x)))
acos(x) = atan2 (sqrt ((1.0 + x) * (1.0 - x)), x)

Es stellt sich heraus, dass die Suche nach einem sinnvollen Beschreibung des CORDIC-iteration für die atan2() für die double ist schwieriger, als ich dachte. Die folgende Webseite wird angezeigt, um eine hinreichend detaillierte Beschreibung, und auch beschreibt zwei alternative Ansätze, die polynomiale approximation und lookup-Tabellen:

http://ch.mathworks.com/examples/matlab-fixed-point-designer/615-calculate-fixed-point-arctangent

Sollte es einfach sein zu addapt Sie den folgenden code, um festen Punkt. Es beschäftigt rationale approximation zur Berechnung des Arkustangens normiert auf den [0 1) - Intervall (können Sie multiplizieren Sie diese mit Pi/2 zu bekommen, die real Arkustangens). Anschließend können Sie bekannte Identitäten zu Holen Sie sich die arcsin/arccos aus dem Arkustangens.

normalized_atan(x) ~ (b x + x^2) /(1 + 2 b x + x^2)

where b = 0.596227

Sich der maximale Fehler ist 0.1620 º

#include <stdint.h>
#include <math.h>

//Approximates atan(x) normalized to the [-1,1] range
//with a maximum error of 0.1620 degrees.

float norm_atan( float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    //Extract the sign bit
    uint32_t ux_s  = sign_mask & (uint32_t &)x;

    //Calculate the arctangent in the first quadrant
    float bx_a = ::fabs( b * x );
    float num = bx_a + x * x;
    float atan_1q = num /( 1.f + bx_a + num );

    //Restore the sign bit
    uint32_t atan_2q = ux_s | (uint32_t &)atan_1q;
    return (float &)atan_2q;
}

//Approximates atan2(y, x) normalized to the [0,4) range
//with a maximum error of 0.1620 degrees

float norm_atan2( float y, float x )
{
    static const uint32_t sign_mask = 0x80000000;
    static const float b = 0.596227f;

    //Extract the sign bits
    uint32_t ux_s  = sign_mask & (uint32_t &)x;
    uint32_t uy_s  = sign_mask & (uint32_t &)y;

    //Determine the quadrant offset
    float q = (float)( ( ~ux_s & uy_s ) >> 29 | ux_s >> 30 ); 

    //Calculate the arctangent in the first quadrant
    float bxy_a = ::fabs( b * x * y );
    float num = bxy_a + y * y;
    float atan_1q =  num /( x * x + bxy_a + num );

    //Translate it to the proper quadrant
    uint32_t uatan_2q = (ux_s ^ uy_s) | (uint32_t &)atan_1q;
    return q + (float &)uatan_2q;
} 

In Fall müssen Sie mehr Präzision, gibt es eine rationale Funktion 3. Ordnung:

normalized_atan(x) ~ ( c x + x^2 + x^3) /( 1 + (c + 1) x + (c + 1) x^2 + x^3)

where c = (1 + sqrt(17)) /8

die eine maximale Annäherung Fehler 0.00811 º

http://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Expression_as_definite_integrals

Könnte man tun, dass die integration numerisch mit Ihrem Wurzel-Funktion, die approximation mit einer unendlichen Reihe:

Verwenden Sie eine Polynom-approximation. Least-squares-fit ist am einfachsten (Microsoft Excel) und Tschebyscheff-approximation genauer ist.

Diese Frage wurde abgedeckt, bevor: Wie Trigonometrische Funktionen funktionieren?