Comment puis-je tenir compte des erreurs d'arrondi dans l'arithmétique en virgule flottante pour les fonctions trigonométriques inverse (et sqrt) (en C)?

Question

I ai une fonction relativement complexe qui prend plusieurs valeurs doubles qui représentent deux vecteurs dans un espace à 3 de la forme (amplitude, latitude, longitude) où la latitude et la longitude sont en radians, et un angle. Le but de cette fonction est de faire tourner le premier vecteur vers le second par l'angle spécifié et renvoyer le vecteur résultant. Je l'ai déjà vérifié que le code est logiquement correct et fonctionne.

Le but attendu de la fonction est pour le graphisme, donc double précision est pas nécessaire; Cependant, sur la plate-forme cible, triglycéride (et sqrt) fonctions qui prennent des flotteurs (sinf, FPVO, atan2f, asinf, acosf et sqrtf spécifiquement) un travail plus rapide sur double que sur flotteurs (probablement parce que l'instruction de calcul de ces valeurs peut effectivement nécessiter une double, si un flotteur est passé, la valeur doit être coulé à un double, ce qui nécessite la copie vers une zone avec plus de mémoire - les frais généraux par exemple). En conséquence, toutes les variables impliquées dans la fonction sont double précision.

Voici la question: Je suis en train d'optimiser ma fonction afin qu'il puisse être appelé plusieurs fois par seconde. J'ai donc remplacé les appels à sin, cos, sqrt, et ainsi de suite avec les appels vers les versions à virgule flottante de ces fonctions, car elles entraînent une augmentation de la vitesse 3-4 fois ensemble. Cela fonctionne pour presque toutes les entrées; Cependant, si des vecteurs d'entrée sont proches parallèles aux vecteurs unitaires type (i, j, k), des erreurs d'arrondi pour les différentes fonctions accumuler suffisamment pour provoquer des appels successifs à sqrtf ou inverse des fonctions trigonométriques (asinf, acosf, atan2f) pour passer des arguments qui sont à peine en dehors du domaine de ces fonctions.

Alors, je reste avec ce dilemme: ou bien je ne peux appeler des fonctions doubles de précision et d'éviter le problème (et finir avec une limite d'environ 1.300.000 opérations vectorielles par seconde), ou je peux essayer de trouver autre chose . En fin de compte, je voudrais un moyen d'assainir l'entrée à l'inverse des fonctions trigonométriques pour prendre en charge des cas de bord (il est trivial pour le faire pour sqrt: abs juste utilisation). Branching est pas une option, car même une seule instruction conditionnelle ajoute tellement frais généraux que les gains de performance sont perdus.

Alors, des idées?

Edit: quelqu'un a exprimé la confusion sur mon double par rapport à l'aide d'opérations en virgule flottante. La fonction est beaucoup plus rapide si je stocke en fait toutes mes valeurs dans des conteneurs à double dimension (par exemple des variables de type double) que si je les entreposer dans des conteneurs float taille. Cependant, les opérations flottantes trigonométriques de précision des points sont plus rapides que la double précision des opérations trigonométriques pour des raisons évidentes.

Answer 1

En gros, vous devez trouver un algorithme numériquement stable qui permet de résoudre votre problème. Il n'y a pas de solutions génériques à ce genre de chose, il doit être fait pour votre cas spécifique en utilisant des concepts tels que le condition numéro si les étapes individuelles. Et il peut en effet être impossible si le problème sous-jacent lui-même est mal conditionné.

Answer 2

simple précision en virgule flottante introduit en soi une erreur. Donc, vous devez construire vos calculs afin que toutes les comparaisons ont un certain degré de « slops » en utilisant un facteur epsilon, et vous avez besoin d'entrées Désinfectez aux fonctions avec des domaines limités.

L'ancien est assez facile lorsque la ramification, par exemple

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

mais c'est en désordre. Blocage des entrées de domaine est un peu plus délicat, mais mieux. La clé est d'utiliser conditionnelle opérateurs de déplacer , qui quelque chose de général do comme

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

dans un seul op, sans encourir une branche.

Ceux-ci varient en fonction de l'architecture. Sur l'unité de virgule flottante x87, vous pouvez le faire avec le FCMOV conditionnel mouvement op , mais c'est maladroite, car elle dépend de drapeaux de condition étant précédemment définie, il est donc lent. En outre, il n'y a pas un compilateur compatible intrinsèque pour cmov. Ceci est l'une des raisons pour lesquelles nous évitons x87 à virgule flottante en faveur de mathématiques de SSE2 si possible.

est conditionnel mouvement beaucoup mieux pris en charge dans l'ESS en associant un opérateur comparaison avec un niveau du bit. Ceci est préférable, même pour les mathématiques scalaire:

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

sont Compilateurs mieux, mais pas génial, à environ émettre ce genre de modèle pour ternaires lorsque les mathématiques scalaire SSE2 est activé. Vous pouvez le faire avec le drapeau du compilateur /arch:sse2 sur MSVC ou -mfpmath=sse sur GCC.

Sur le PowerPC et bien d'autres architectures RISC, fsel() est un opcode matériel et donc généralement un intrinsèque du compilateur ainsi.

Answer 3

Avez-vous regardé Graphics Programmation Livre noir ou la remise peut-être les calculs hors de votre GPU?