Come posso spiegare errori di arrotondamento in virgola mobile per la trigonometria inversa (e sqrt) funzioni (in C)?

Question

Ho una funzione piuttosto complicato che richiede diversi valori doppi che rappresentano due vettori in 3-spazio del modulo (ampiezza, latitudine, longitudine) dove latitudine e longitudine sono in radianti, e un angolo. Lo scopo della funzione è ruotare il primo vettore attorno al secondo un angolo specificato e restituire il vettore risultante. Ho già verificato che il codice è logicamente corretto e lavora.

Lo scopo previsto della funzione è per la grafica, quindi doppia precisione non è necessaria; tuttavia, sulla piattaforma di destinazione, trig (e sqrt) le funzioni che accettano i galleggianti (sinf, cosf, atan2f, asinf, acosf e sqrtf in particolare) il lavoro più veloce su doppie che su carri (probabilmente perché le istruzioni per calcolare tali valori può effettivamente richiedere una doppio, se è passato un galleggiante, il valore deve essere fuso ad una doppia, che richiede copiandolo in una zona con più memoria - ossia dall'alto). Di conseguenza, tutte le variabili coinvolte nella funzione sono doppia precisione.

Ecco il problema: Sto cercando di ottimizzare la mia funzione in modo che possa essere chiamato più volte al secondo. Ho pertanto sostituito le chiamate a sin, cos, sqrt, eccetera di chiamate verso le versioni virgola mobile di tali funzioni, quali risultano in un aumento di velocità 3-4 volte complessiva. Questo funziona per quasi tutti gli ingressi; tuttavia, se i vettori di ingresso sono vicini a parallelo con i versori normale (i, j, o k), arrotondamento errori per le varie funzioni costruire sufficiente a causare successivamente chiamate a sqrtf o inversa funzioni trigonometriche (asinf, acosf, atan2f) per passare gli argomenti che sono appena al di fuori del dominio di tali funzioni.

Quindi, io sono rimasto con questo dilemma: o io posso chiamare solo funzioni doppie di precisione ed evitare il problema (e finire con un limite di circa 1.300.000 operazioni vettoriali al secondo), o posso provare a venire con qualcos'altro . In definitiva, mi piacerebbe un modo per sterilizzare l'input per l'inverso funzioni trigonometriche per prendersi cura di casi limite (è banale per farlo per sqrt: solo uso abs). Branching non è un'opzione, come anche una sola dichiarazione condizionale aggiunge così tanto in testa che eventuali miglioramenti delle prestazioni sono persi.

Quindi, tutte le idee?

Edit: qualcuno ha espresso la confusione sopra i miei utilizzando doppie rispetto a operazioni in virgola mobile. La funzione è molto più veloce se ho effettivamente memorizzare tutti i miei valori in contenitori di dimensioni doppie (variabili OSSIA doppio tipo) che se li devo conservare in contenitori float-size. Tuttavia, operazioni in virgola mobile di precisione trigonometriche sono più veloci di doppia precisione operazioni trigonometriche per ovvi motivi.

Answer 1

In sostanza, è necessario trovare un href="http://en.wikipedia.org/wiki/Numerical_stability" rel="nofollow"> algoritmo di numericamente stabile numero di condizione se le singole fasi. E può infatti essere impossibile se il problema alla base stessa è mal condizionata.

Answer 2

virgola mobile a singola precisione introduce intrinsecamente errore. Quindi, è necessario costruire la matematica in modo che tutti i confronti hanno un certo grado di "slop" utilizzando un fattore Epsilon, e avete bisogno di ingressi Disinfettare a funzioni con i domini limitate.

Il primo è abbastanza facile quando ramificazione, ad esempio

bool IsAlmostEqual( float a, float b ) { return fabs(a-b) < 0.001f; } // or
bool IsAlmostEqual( float a, float b ) { return fabs(a-b) < (a * 0.0001f); } // for relative error

Ma questo è disordinato. Serraggio ingressi dominio è un po 'più complicato, ma meglio. La chiave è quella di utilizzare condizionale operatori spostare , che in do qualcosa di generale come

float ExampleOfConditionalMoveIntrinsic( float comparand, float a, float b ) 
{ return comparand >= 0.0f ? a : b ; }

in un unico op, senza incorrere in un ramo.

Questi variano a seconda dell'architettura. Sull'unità virgola mobile x87 si può fare con la FCMOV condizionale-move op , ma che è maldestra perché dipende flag di condizione di essere impostato in precedenza, quindi è lento. Inoltre, non c'è un intrinseco compilatore coerente per cmov. Questo è uno dei motivi per cui evitiamo x87 floating point a favore di SSE2 scalare la matematica, se possibile.

mossa Condizionale sarà molto meglio supportato in SSE di abbinare un confronto operatore con un AND bit per bit. Questo è ancora preferibile per la matematica scalare:

// assuming you've already used _mm_load_ss to load your floats onto registers 
__m128 fsel( __m128 comparand, __m128 a, __m128 b ) 
{
    __m128 zero = {0,0,0,0};
    // set low word of mask to all 1s if comparand > 0
    __m128 mask = _mm_cmpgt_ss( comparand, zero );  
    a = _mm_and_ss( a, mask );    // a = a & mask 
    b = _mm_andnot_ss( mask, b ); // b = ~mask & b
    return _mm_or_ss( a, b );     // return a | b
    }
}

I compilatori sono meglio, ma non ottima, a circa emettono questo tipo di modello per ternari quando SSE2 scalare la matematica è abilitato. Potete farlo con il compilatore /arch:sse2 bandiera sul MSVC o -mfpmath=sse su GCC.

Sulla PowerPC e molte altre architetture RISC, fsel() è un opcode hardware e quindi di solito un'intrinseca compilatore pure.

Answer 3

Hai guardato la grafica di programmazione Black Book o forse consegnando i calcoli fuori alla vostra GPU?