Gibt es eine genaue Annäherung der acos() - Funktion?

Question

Ich brauche eine acos() Funktion mit doppelter Genauigkeit innerhalb eines Compute-Shaders. Da es keine eingebaute Funktion von acos() in GLSL mit doppelter Genauigkeit gibt, habe ich versucht, meine eigene zu implementieren.

Zuerst habe ich eine Taylor-Serie wie die Gleichung aus Wiki - Taylor series mit vorberechneten Fakultät Werte implementiert. Aber das scheint um 1 herum ungenau zu sein. Der maximale Fehler betrug ungefähr 0,08 mit 40 Iterationen.

Ich implementiert auch this method, die sehr gut auf CPU mit einem maximalen Fehler von -2.22045e-16 funktioniert, aber ich habe einige Probleme, dies innerhalb des Shaders zu implementieren.

Derzeit verwende ich eine acos() Approximationsfunktion von here, wo jemand seine Approximationsfunktionen auf this Seite veröffentlicht. Ich verwende die genaueste Funktion dieser Seite und jetzt bekomme ich einen maximalen Fehler von -7.60454e-08, aber auch dieser Fehler ist zu hoch.

Mein Code dieser Funktion ist:

double myACOS(double x) 
{ 
    double part[4]; 
    part[0] = 32768.0/2835.0*sqrt(2.0-sqrt(2.0+sqrt(2.0+sqrt(2.0+2.0*x)))); 
    part[1] = 256.0/135.0*sqrt(2.0-sqrt(2.0+sqrt(2.0+2.0*x))); 
    part[2] = 8.0/135.0*sqrt(2.0-sqrt(2.0+2.0*x)); 
    part[3] = 1.0/2835.0*sqrt(2.0-2.0*x); 
    return (part[0]-part[1]+part[2]-part[3]); 
} 

Kennt jemand eine andere Implementierungsmethode von acos(), die sehr genau ist und -wenn möglich- einfach in einem Shader zu implementieren?

Einige Systeminformationen:

• Nvidia GT 555M
• laufen OpenGL 4.3 mit optirun

Answer 1

Meine aktuellen genauen Shader-Implementierung von 'acos()' ist eine Mischung aus dem üblichen Taylor aus Serie und die Antwort von Bence. Mit 40 Iterationen erhalte ich eine Genauigkeit von 4.44089e-16 gegenüber der 'acos()' Implementierung von math.h. Vielleicht ist es nicht die beste, aber es funktioniert für mich:

Und hier ist es:

double myASIN2(double x) 
{ 
    double sum, tempExp; 
    tempExp = x; 
    double factor = 1.0; 
    double divisor = 1.0; 
    sum = x; 
    for(int i = 0; i < 40; i++) 
    { 
     tempExp *= x*x; 
     divisor += 2.0; 
     factor *= (2.0*double(i) + 1.0)/((double(i)+1.0)*2.0); 
     sum += factor*tempExp/divisor; 
    } 
    return sum; 
} 

double myASIN(double x) 
{ 
    if(abs(x) <= 0.71) 
     return myASIN2(x); 
    else if(x > 0) 
     return (PI/2.0-myASIN2(sqrt(1.0-(x*x)))); 
    else //x < 0 or x is NaN 
     return (myASIN2(sqrt(1.0-(x*x)))-PI/2.0); 

} 

double myACOS(double x) 
{ 
    return (PI/2.0 - myASIN(x)); 
} 

Kommentare, was besser gemacht werden könnte? Zum Beispiel mit einer LUT für die Werte von Faktor, aber in meinem Shader 'acos()' wird nur einmal aufgerufen, so dass es nicht nötig ist.

Answer 2

Die NVIDIA GT 555M GPU ist ein Gerät mit Rechenkapazität 2.1, daher gibt es native Hardwareunterstützung für grundlegende Operationen mit doppelter Genauigkeit, einschließlich fused multipy-add (FMA). Wie bei allen NVIDIA-GPUs wird die Quadratwurzeloperation emuliert. Ich kenne CUDA, aber nicht GLSL. Gemäß der Version 4.3 der GLSL specification legt sie FMA mit doppelter Genauigkeit als Funktion fma() offen und stellt eine Quadratwurzel mit doppelter Genauigkeit, sqrt(), bereit. Es ist nicht klar, ob die sqrt() Implementierung korrekt nach IEEE-754 Regeln gerundet wird. Ich gehe davon aus, dass es in Analogie zu CUDA ist.

Anstatt eine Taylor-Serie zu verwenden, würde man ein Polynom verwenden und damit die Anzahl der benötigten Terme reduzieren. Minimax-Approximationen werden typischerweise unter Verwendung einer Variante der Remez algorithm erzeugt. Um Geschwindigkeit und Genauigkeit zu optimieren, ist der Einsatz von FMA unerlässlich. Die Auswertung des Polynoms mit der Horner scheme ist günstig für eine hohe Anwachsung. In dem folgenden Code wird ein Horner-Schema zweiter Ordnung verwendet.Wie in DanceIgels answer wird acos bequem unter Verwendung einer asin Approximation als der Grundbaustein in Verbindung mit mathematischen Standardidentitäten berechnet.

Bei 400M-Testvektoren war der maximale relative Fehler, der mit dem unten stehenden Code zu sehen war, 2,67e-16, während der maximal beobachtete Fehler ulp 1,442 ul ist.

/* compute arcsin (a) for a in [-9/16, 9/16] */ 
double asin_core (double a) 
{ 
    double q, r, s, t; 

    s = a * a; 
    q = s * s; 
    r =    5.5579749017470502e-2; 
    t =   -6.2027913464120114e-2; 
    r = fma (r, q, 5.4224464349245036e-2); 
    t = fma (t, q, -1.1326992890324464e-2); 
    r = fma (r, q, 1.5268872539397656e-2); 
    t = fma (t, q, 1.0493798473372081e-2); 
    r = fma (r, q, 1.4106045900607047e-2); 
    t = fma (t, q, 1.7339776384962050e-2); 
    r = fma (r, q, 2.2372961589651054e-2); 
    t = fma (t, q, 3.0381912707941005e-2); 
    r = fma (r, q, 4.4642857881094775e-2); 
    t = fma (t, q, 7.4999999991367292e-2); 
    r = fma (r, s, t); 
    r = fma (r, s, 1.6666666666670193e-1); 
    t = a * s; 
    r = fma (r, t, a); 

    return r; 
} 

/* Compute arccosine (a), maximum error observed: 1.4316 ulp 
    Double-precision factorization of π courtesy of Tor Myklebust 
*/ 
double my_acos (double a) 
{ 
    double r; 

    r = (a > 0.0) ? -a : a; // avoid modifying the "sign" of NaNs 
    if (r > -0.5625) { 
     /* arccos(x) = pi/2 - arcsin(x) */ 
     r = fma (9.3282184640716537e-1, 1.6839188885261840e+0, asin_core (r)); 
    } else { 
     /* arccos(x) = 2 * arcsin (sqrt ((1-x)/2)) */ 
     r = 2.0 * asin_core (sqrt (fma (0.5, r, 0.5))); 
    } 
    if (!(a > 0.0) && (a >= -1.0)) { // avoid modifying the "sign" of NaNs 
     /* arccos (-x) = pi - arccos(x) */ 
     r = fma (1.8656436928143307e+0, 1.6839188885261840e+0, -r); 
    } 
    return r; 
}