Racine inverse - existe en C++ ?

Version imprimable

23/12/2009, 20h33
Nanoc

Racine inverse - existe en C++ ?

Bonsoir,

J'ai un programme qui doit calculer entre autre la racine carrée inverse d'un nombre en simple ou double précision. Soit:

Code:

1. / sqrt(x);

et ceci des millions de fois. Des mesures avec callgrind montrent que le programme passe en effet près de 25% de son temps dans cette opération.

La question que je me pose est de savoir si il existe une fonction permettant d'obtenir directement 1/sqrt(x) sans passer par sqrt(x) et une division ?
Ou alors est-ce que le compilateur est capable d'optimiser cela en utilisant des instructions processeurs ?

Auriez-vous des informations à ce sujet ?

Merci d'avance !

P.S.: Je ne cherche PAS des trucs approchés comme on trouve dans les jeux vidéo. J'ai besoin d'une solution aussi exacte que 1./sqrt(x).
23/12/2009, 21h45
Davidbrcz

Ce PDF devrait de brancher. Attention, ce ne sont pas des maths de bas niveau.

Citation:

P.S.: Je ne cherche PAS des trucs approchés comme on trouve dans les jeux vidéo. J'ai besoin d'une solution aussi exacte que 1./sqrt(x).

Va faire des maths alors :D. En informatique, tu auras forcément une valeur approchée. Peut être à 10^(-10), 10^(-100) mais elle reste approchée.
24/12/2009, 01h18
deadalnix

La solution n'est précise qu'a 10^-3 près environs. Ce n'est donc pas ce que cherche notre homme.

A priori, il est existe des instructions processeur pour faire cela. Pourquoi ne regardes-tu pas le code généré par le compilo afin de savoir si tu ne peux pas résoudre ton problème par des réglages compilo ? (pour g++ -O3 et -ffast-math)
24/12/2009, 10h33
Davidbrcz

Citation:

Envoyé par deadalnix

La solution n'est précise qu'a 10^-3 près environs. Ce n'est donc pas ce que cherche notre homme.

A priori, il est existe des instructions processeur pour faire cela. Pourquoi ne regardes-tu pas le code généré par le compilo afin de savoir si tu ne peux pas résoudre ton problème par des réglages compilo ? (pour g++ -O3 et -ffast-math)

La méthode proposée se base sur l'extraction de racines par la méthode de Newton, qui est un algorithme à convergence quadratique, donc si la précision n'est pas suffisante, il suffit de recommencer :D
03/01/2010, 04h30
Mac LAK

Question bête : sachant que 1/sqrt(x) est égal à (x^(-1/2)), as-tu essayé un appel à pow avec -0.5 en guise d'exposant, pour voir ce que ça donne en terme de performances ?

Sinon, faudrait peut-être regarder du côté des instructions SSE. Mais de mémoire, il n'y a pas d'opération native du FPU qui effectue directement 1/sqrt(x).
03/01/2010, 11h22
screetch

_mm_rsqrt_ps en SSE
http://msdn.microsoft.com/en-us/libr...8VS.71%29.aspx
sous GCC tu pourrais peut etre augmenter les perfs avec les options
'-mfpmath=sse', '-msse2'
03/01/2010, 11h27
loufoque

Citation:

sous GCC tu pourrais peut etre augmenter les perfs avec les options
'-mfpmath=sse', '-msse2'

Il faut mieux donner le bon -march et laisser le compilo choisir le mieux
03/01/2010, 11h31
mohira

Bonjour,
Vous pourriez peut-être essayer de programmer ceci pour un calcul de 1/sqrt(x): si on considère que y est une valeur approchée,on pose (y+a)^2=1./x et on approxime par y^2+2*a*y=1./x d'où a = (1./(x*y)-y)/2.,on fait y=y+a,et on itère avec un test de convergence.
A la question de l'initialisation de y,je pense que,si on sait identifier l'exposant et la mantisse de x,on pourrait prendre la même mantisse pour y(initial) et pour exposant l'exposant de x divisé par -2.
03/01/2010, 15h01
deadalnix

Citation:

Envoyé par Davidbrcz

La méthode proposée se base sur l'extraction de racines par la méthode de Newton, qui est un algorithme à convergence quadratique, donc si la précision n'est pas suffisante, il suffit de recommencer :D

Ça sera de toute façon moins rapide et moins précis que l'instruction cablée.
03/01/2010, 23h08
mohira

Il est certain que,si l'instruction est câblée,la discussion sur une méthode logicielle n'a plus lieu.Bien sûr,quand je pensais à une méthode logicielle,c'est parce que je pensais qu'on supposait a priori que ce n'était pas le cas.C'était l'intérêt de la question initiale.Si la fonction sqrt() est câblée,la discussion est close d'elle-même.
04/01/2010, 12h37
Lavock

La question serait plutôt amha pourquoi tu passe 25% de ton temps dans 4 lignes d'assembleur ?
08/01/2010, 21h42
RaphaelSP

Citation:

Envoyé par screetch

_mm_rsqrt_ps en SSE
http://msdn.microsoft.com/en-us/libr...8VS.71%29.aspx
sous GCC tu pourrais peut etre augmenter les perfs avec les options
'-mfpmath=sse', '-msse2'

Cette instruction fait en effet le travail, mais elle est relativement imprécise: environ 11 bits sur 24 sont significatifs.
Il est donc nécessaire soit de refaire une itération de Newton-Raphson, soit de passer par la racine plus précise puis par la division (_mm_sqrt_ps puis _mm_div_ps) (ce qui est beaucoup plus long - cf. guide d'optimisation Intel: http://www.intel.com/products/proces...uals/index.htm).
10/01/2010, 11h40
Nanoc

Merci pour vos réponses.

pow(x,-0.5); est plus lent et de beaucoup.

Pour les instructions SSE, je n'ai pas vu de différence notable. Je ne sais pas si le compilateur ne les activait pas par défaut.

Citation:

Envoyé par Nanoc

Pour les instructions SSE, je n'ai pas vu de différence notable. Je ne sais pas si le compilateur ne les activait pas par défaut.

C'est bizarre, non ?
Je n'y connais pas grand chose en SSE, mais a priori ces instructions permettent de faire quatre racine inverse en une passe.

Avec le test naïf ci-dessous, avec vs2010, j'obtiens bien un boost x4 :

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
 
#include <windows.h>
#include <iostream>
 
#include <cmath>
#include "xmmintrin.h"
 
#define LENGTH 100000000
 
void TestInvSqrt()
{
   double d = 0;
   for(int i = 1 ; i < LENGTH ; i++)
   {
      d += 1/sqrt((double)i);
   }
   printf("%f\n", d);
}
 
 
void TestInvSqrt2()
{
   double d = 0;
   for(int i = 1 ; i < LENGTH ; i+=4)
   {
      __declspec(align(16)) float tab[4] = {i, i +1, i+2, i+3};
      __m128* m128 = (__m128*)tab;
      *m128 = _mm_rsqrt_ps(*m128);
 
      d += tab[0] + tab[1] + tab[2] + tab[3];
   }
   printf("%f\n", d);
}
 
 
int main()
{
   int time  = GetTickCount();
   TestInvSqrt();
   std::cout << GetTickCount() - time << std::endl; 
 
   time  = GetTickCount();
   TestInvSqrt2();
   std::cout << GetTickCount() - time << std::endl;
}

11/01/2010, 14h14
RaphaelSP
Citation:
Envoyé par Arzar

C'est bizarre, non ?
Je n'y connais pas grand chose en SSE, mais a priori ces instructions permettent de faire quatre racine inverse en une passe.

Avec le test naïf ci-dessous, avec vs2010, j'obtiens bien un boost x4 :

Code:

1 2 //...
Nanoc, comme vient de le montrer Arzar, c'est en utilisant explicitement les instructions SIMD (version packed, pas scalaire) que tu obtiendras vraiment des gains.

Activer les intrinsèques dans les options permet au compilateur de les utiliser, mais il ne parallélisera pas pour toi - à moins d'utiliser ICC (un de leurs arguments de vente est la parallélisation automatique).

edit: Note: comme je le précisais plus haut, utiliser seulement _mm_rsqrt_ps ne te permet d'obtenir que 12 bits de précision. Si ça te suffit, tant mieux. Sinon, un raffinage reste à faire.

edit2: mais de toute façon, si tu passes en SIMD, tu peux probablement commencer avec les instructions câblées et précises. Si le reste de ton calcul est parallélisable aussi, tu as des gains bien plus grand à espérer.
11/01/2010, 17h30
Lavock

Je vais finir par croire qu'ils ont raison, un peu de théorie ne fait pas de mal >< !

SSE, historiquement, c'est pour de la 3D. Du coup, cela PERMET le parallélisme, ou pour être plus exact, de la vectorisation. A ma connaissance, GCC depuis GCC 4 introduit la vectorisation automatique pour SSE et SSE2 >< !

Après, si c'est vectorisable, en fonction de comment, soit on fait les opérations sur GPU (OpenCL), soit on se contente des instructions SSE... Mais comme déjà dit, GCC le fait pour toi normalement. Visual, je sais pas.

Reste la façon d'activer les instruction SSE. C'est un truc qui m'ennuie au plus au point personnellement.

Si tu utilise un -march approprié, tu les actives (car elles sous-entendent -msse[x]).

Tu peux aussi utilise -msse[x], histoire de sélectionné non pas une architecture, mais juste les truc requis.

La doc dis qu'il faut mettre -mpfmath=sse pour les utiliser automatiquement, mais sans ça, il les utilise quand même >< ??? (c'est ça le truc ennuyant car pas clair du tout).

Une chose dont je suis pas certain, c'est que la version "imprécise" nécessite -funsafe-math-optimizations. Je dirais que c'est ce qui me parait logique.

Bon ça résout pas tout. L'important étant surtout de savoir si c'est normal que tu fasse autant de division.
11/01/2010, 18h09
Mac LAK

Citation:

Envoyé par Lavock

Visual, je sais pas.

Visual utilise surtout les SSE (de façon automatique) pour optimiser les calculs en flottant et/ou entiers, mais pas spécialement pour paralléliser à ta place : http://msdn.microsoft.com/fr-fr/library/7t5yh4fd.aspx

Pour le reste, VS fournit les fonctions intrinsèques de commande directe des copros SSE, ce qui permet de construire directement des boucles parallèles utilisant ces registres.