Esiste un'accurata approssimazione della funzione acos()?

Question

Ho bisogno di una funzione acos() con una precisione doppia all'interno di uno shader di elaborazione. Poiché non esiste una funzione integrata di acos() in GLSL con una precisione doppia, ho provato a implementare la mia.

Inizialmente, ho implementato una serie di Taylor come l'equazione da Wiki - Taylor series con valori di facoltà precalcolati. Ma questo sembra essere inaccurato intorno a 1. L'errore massimo era qualcosa di circa 0.08 con 40 iterazioni.

Ho anche implementato this method che funziona molto bene sulla CPU con un errore massimo di -2.22045e-16, ma ho alcuni problemi per implementarlo nello shader.

Attualmente, sto usando una funzione di approssimazione acos() da here dove qualcuno ha pubblicato le sue funzioni di approssimazione sul sito this. Sto usando la funzione più precisa di questo sito e ora ottengo un errore massimo di -7.60454e-08, ma anche che l'errore è troppo alto.

Il mio codice di questa funzione è:

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]); 
} 

Qualcuno sa un altro metodo attuazione acos() che è molto accurata e -se possibile- facile da implementare in uno shader?

alcune informazioni di sistema:

• Nvidia GT 555M
• esecuzione OpenGL 4.3 con optirun

Answer 1

mio attuale implementazione di shader accurata di 'acos()' è un mix fuori il solito Taylor serie e la risposta da Bence. Con 40 iterazioni ottengo una precisione di 4.44089e-16 per l'implementazione 'acos()' da math.h. Forse non è la migliore, ma funziona per me:

E qui è:

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)); 
} 

Eventuali commenti, cosa si poteva fare di meglio? Ad esempio, usando una LUT per i valori di factor, ma nel mio shader 'acos()' viene chiamato una sola volta, quindi non ce n'è bisogno.

Answer 2

La GPU NVIDIA GT 555M è un dispositivo con capacità di calcolo 2.1, pertanto è disponibile il supporto hardware nativo per operazioni di precisione doppia di base, tra cui fused multipy-add (FMA). Come con tutte le GPU NVIDIA, l'operazione di radice quadrata viene emulata. Conosco CUDA, ma non GLSL. In base alla versione 4.3 di GLSL specification, espone FMA a doppia precisione come funzione fma() e fornisce una radice quadrata a precisione doppia, sqrt(). Non è chiaro se l'implementazione sqrt() sia arrotondata correttamente in base alle regole IEEE-754. Presumo che lo sia, per analogia con CUDA.

Invece di utilizzare una serie di Taylor, si vorrebbe usare un polinomio minimax approximation, riducendo così il numero di termini richiesti. Le approssimazioni Minimax vengono generalmente generate utilizzando una variante di Remez algorithm. Per ottimizzare sia la velocità che la precisione, l'uso di FMA è essenziale. La valutazione del polinomio con Horner scheme è conduttiva per l'alta accrastazione. Nel codice qui sotto, viene utilizzato uno schema Horner di secondo ordine.Come in DanceIgel answer, lo acos viene calcolato convenientemente utilizzando un'approssimazione asin come base elementare in combinazione con le identità matematiche standard.

Con i vettori di test da 400 M, l'errore relativo massimo visto con il codice sottostante era 2,67e-16, mentre l'errore massimo di ulp osservato è 1,444 ulp.

/* 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; 
}