Comment normaliser 4 vecteurs rapidement via les instructions SIMD

Une grande partie des performances dégagées par notre moteur sont dues au concept de SIMD; Simple Instruction Multiple Data. Ce sont des instructions processeurs (SSE, SSE2, SSE3, SSE4a, 3dNow!) qui permettent d'accélérer des calculs bas niveau.


L'exemple suivant montre comment faire une addition de vecteur classique

 struct vector  
 {  
   float x,y,z;  
 };  
 vector add(const vector &a, const vector &b)  
 {  
   vector c;  
   c.x = a.x + b.x;  
   c.y = a.y + b.y;  
   c.z = a.z + b.z;  
   return c;  
 }  

Et maintenant la version SIMD:

 # include <xmmintrin.h>  
 struct vector  
 {  
   __m128 m;  
   // Pour acceder aux composantes x, y, z  
   // vous pouvez simplement faire un cast ((float*)&m)[ 0 1 2 ou 3]  
 };  
 vector add(const vector &a, const vector &b)  
 {  
   vector c;  
   c.m = _mm_add_ps(a.m, b.m);  
   return c;  
 }  

Entre les deux versions, le facteur de speedup se situe entre 3 et 4. Si vous développez des applications nécessitants de nombreux calculs 3D, utilisez absolument ces instructions!


Voici maintenant une astuce pour gagner du temps sur les normalisations.


Il est plus performant de multiplier par la sqrt inverse (_mm_rsqrt_ps) que de diviser par la sqrt. En revanche la fonction _mm_rsqrt_ps a une précision moindre. Pour combler ce manque de précision, on fait une nouvelle itération dans le calcul de la racine, appelée "Iteration de Newton-Raphson".

 __m128    fastrsqrt( const __m128 v )  
 {  
       static const __m128 _half4 = _mm_set_ps1( 0.5f );  
       static const __m128 _three = _mm_set_ps1( 3.0f );  
       const __m128 nr = _mm_rsqrt_ps( v ); // Newton-Raphson iteration  
       const __m128 muls = _mm_mul_ps( _mm_mul_ps( v, nr ), nr );  
       return _mm_mul_ps( _mm_mul_ps( _half4, nr ), _mm_sub_ps( _three, muls ) );  
 }  
 void Normalize4(vector p[4])  
 {  
     // Store x, y, and z separately in 3 __m128  
     __m128 x4 = _mm_set_ps(p[3][0], p[2][0], p[1][0], p[0][0]);  
     __m128 y4 = _mm_set_ps(p[3][1], p[2][1], p[1][1], p[0][1]);  
     __m128 z4 = _mm_set_ps(p[3][2], p[2][2], p[1][2], p[0][2]);  
     x4 = _mm_mul_ps(x4, x4); // x * x  
     y4 = _mm_mul_ps(y4, y4); // y * y  
     z4 = _mm_mul_ps(z4, z4); // z * z  
     const vector l = fastrsqrt(_mm_add_ps(z4, _mm_add_ps(y4, x4))); // sqrt(x*x+y*y+z*z)  
     p[0].m = _mm_mul_ps(p[0].m, _mm_set_ps1(l[0])); // load results in packet  
     p[1].m = _mm_mul_ps(p[1].m, _mm_set_ps1(l[1]));  
     p[2].m = _mm_mul_ps(p[2].m, _mm_set_ps1(l[2]));  
     p[3].m = _mm_mul_ps(p[3].m, _mm_set_ps1(l[3]));  
 }  

Si vous avez aimé l'article n'hésitez pas à laisser un commentaire, ou me suivre sur twitter

Immersion Engine dans la presse!

Nous sommes dans le magazine Programmez du mois de Janvier!
A la page numéro 24, vous trouverez une page parlant de notre projet.

Je vous invite à acheter le magazine pour découvrir les autres passionnants projets de fin d'études présenté par les étudiants de notre promotion. Vous pouvez vous le procurer ici.

Bonne lecture!

Si vous avez aimé l'article n'hésitez pas à laisser un commentaire, ou me suivre sur twitter