Pourquoi le scalaire SSE sqrt (x) est-il plus lent que rsqrt (x) * x?

106

J'ai profilé certaines de nos mathématiques de base sur un Intel Core Duo, et en examinant diverses approches de la racine carrée, j'ai remarqué quelque chose d'étrange: en utilisant les opérations scalaires SSE, il est plus rapide de prendre une racine carrée réciproque et de la multiplier. pour obtenir le sqrt, que pour utiliser l'opcode sqrt natif!

Je le teste avec une boucle quelque chose comme:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

J'ai essayé cela avec quelques corps différents pour la TestSqrtFunction, et j'ai des timings qui me grattent vraiment la tête. Le pire de tous était de loin d'utiliser la fonction native sqrt () et de laisser le compilateur «intelligent» «optimiser». À 24ns / float, en utilisant le FPU x87, c'était pathétiquement mauvais:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

La prochaine chose que j'ai essayée était d'utiliser un intrinsèque pour forcer le compilateur à utiliser l'opcode scalar sqrt de SSE:

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

C'était mieux, à 11,9 ns / flotteur. J'ai aussi essayé la technique d'approximation farfelue Newton-Raphson de Carmack , qui fonctionnait encore mieux que le matériel, à 4,3 ns / flottant, bien qu'avec une erreur de 1 sur 2 10 (ce qui est trop pour mes besoins).

Le problème, c'était quand j'ai essayé l'opération SSE pour la racine carrée réciproque , puis utilisé une multiplication pour obtenir la racine carrée (x * 1 / √x = √x). Même si cela nécessite deux opérations dépendantes, c'était de loin la solution la plus rapide, à 1,24 ns / flottant et précise à 2-14 :

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Ma question est essentiellement ce qui donne ? Pourquoi l'opcode racine carrée intégré au matériel de SSE est- il plus lent que de le synthétiser à partir de deux autres opérations mathématiques?

Je suis sûr que c'est vraiment le coût de l'opération elle-même, car j'ai vérifié:

  • Toutes les données tiennent dans le cache et les accès sont séquentiels
  • les fonctions sont intégrées
  • dérouler la boucle ne fait aucune différence
  • les indicateurs du compilateur sont définis sur l'optimisation complète (et l'assemblage est bon, j'ai vérifié)

( modifier : la stephentyrone souligne correctement que les opérations sur de longues chaînes de nombres devraient utiliser les opérations de vectorisation SIMD compactées, comme rsqrtps- mais la structure de données du tableau ici est uniquement à des fins de test: ce que j'essaie vraiment de mesurer, ce sont les performances scalaires à utiliser dans le code qui ne peut pas être vectorisé.)

Crashworks
la source
13
x / sqrt (x) = sqrt (x). Ou, en d'autres termes: x ^ 1 * x ^ (- 1/2) = x ^ (1 - 1/2) = x ^ (1/2) = sqrt (x)
Crashworks
6
bien sûr, inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; }. Mais c'est une mauvaise idée car cela peut facilement induire un décrochage du magasin de chargement si le processeur écrit les flottants dans la pile puis les lit immédiatement - jonglant du registre vectoriel à un registre flottant pour la valeur de retour en particulier est une mauvaise nouvelle. En outre, les opcodes machine sous-jacents que les intrinsèques SSE représentent prennent de toute façon des opérandes d'adresse.
Crashworks
4
L'importance de LHS dépend de la génération particulière et du pas d'un x86 donné: mon expérience est que sur tout ce qui va jusqu'à i7, le déplacement des données entre les ensembles de registres (par exemple FPU vers SSE vers eax) est très mauvais, tandis qu'un aller-retour entre xmm0 et pile et le retour ne l'est pas, à cause du transfert de magasin d'Intel. Vous pouvez le chronométrer vous-même pour en être sûr. Généralement, le moyen le plus simple de voir le potentiel LHS est de regarder l'assemblage émis et de voir où les données sont jonglées entre les ensembles de registres; votre compilateur peut faire la chose intelligente, ou pas. En ce qui concerne la normalisation des vecteurs, j'ai écrit mes résultats ici: bit.ly/9W5zoU
Crashworks
2
Pour le PowerPC, oui: IBM dispose d'un simulateur de processeur capable de prédire LHS et de nombreuses autres bulles de pipeline grâce à une analyse statique. Certains PPC ont également un compteur matériel pour LHS que vous pouvez interroger. C'est plus difficile pour le x86; les bons outils de profilage sont plus rares (VTune est quelque peu cassé ces jours-ci) et les pipelines réorganisés sont moins déterministes. Vous pouvez essayer de le mesurer de manière empirique en mesurant les instructions par cycle, ce qui peut être fait avec précision avec les compteurs de performances matérielles. Les registres "instructions retirées" et "cycles totaux" peuvent être lus avec par exemple PAPI ou PerfSuite ( bit.ly/an6cMt ).
Crashworks
2
Vous pouvez également simplement écrire quelques permutations sur une fonction et les chronométrer pour voir si l'une d'entre elles souffre particulièrement de calages. Intel ne publie pas beaucoup de détails sur le fonctionnement de leurs pipelines (qu'ils LHS du tout est en quelque sorte un sale secret), donc une grande partie de ce que j'ai appris a été en regardant un scénario qui provoque un blocage sur d'autres arches (par exemple PPC ), puis en construisant une expérience contrôlée pour voir si le x86 l'a également.
Crashworks

Réponses:

216

sqrtssdonne un résultat correctement arrondi. rsqrtssdonne une approximation de l'inverse, précise à environ 11 bits.

sqrtssgénère un résultat beaucoup plus précis, lorsque la précision est requise. rsqrtssexiste pour les cas où une approximation suffit, mais la rapidité est requise. Si vous lisez la documentation d'Intel, vous trouverez également une séquence d'instructions (approximation réciproque de la racine carrée suivie d'un seul pas de Newton-Raphson) qui donne une précision presque totale (~ 23 bits de précision, si je me souviens bien), et est encore un peu plus rapide que sqrtss.

edit: Si la vitesse est critique, et que vous l'appelez vraiment en boucle pour de nombreuses valeurs, vous devriez utiliser les versions vectorisées de ces instructions, rsqrtpsou les sqrtpsdeux traitent quatre flottants par instruction.

Stephen Canon
la source
3
Le pas n / r vous donne 22 bits de précision (il le double); 23 bits serait exactement une précision totale.
Jasper Bekkers
7
@Jasper Bekkers: Non, ce ne serait pas le cas. Tout d'abord, float a 24 bits de précision. Deuxièmement, sqrtssest correctement arrondi , ce qui nécessite ~ 50 bits avant l'arrondi, et ne peut pas être obtenu en utilisant une simple itération N / R en simple précision.
Stephen Canon
1
C'est certainement la raison. Pour étendre ce résultat: le projet Embree d'Intel ( software.intel.com/en-us/articles/… ), utilise la vectorisation pour ses mathématiques. Vous pouvez télécharger la source sur ce lien et voir comment ils font leurs vecteurs 3/4 D. Leur normalisation vectorielle utilise rsqrt suivi d'une itération de newton-raphson, qui est alors très précis et toujours plus rapide que 1 / ssqrt!
Brandon Pelfrey
7
Une petite mise en garde: x rsqrt (x) donne NaN si x est égal à zéro ou à l'infini. 0 * rsqrt (0) = 0 * INF = NaN. INF rsqrt (INF) = INF * 0 = NaN. Pour cette raison, CUDA sur les GPU NVIDIA calcule des racines carrées à simple précision approximatives sous forme de récip (rsqrt (x)), le matériel fournissant à la fois une approximation rapide de la racine carrée réciproque et réciproque. Évidemment, des contrôles explicites traitant les deux cas particuliers sont également possibles (mais seraient plus lents sur le GPU).
njuffa
@BrandonPelfrey Dans quel fichier avez-vous trouvé l'étape Newton Rhapson?
fredoverflow
7

Ceci est également vrai pour la division. MULSS (a, RCPSS (b)) est bien plus rapide que DIVSS (a, b). En fait, il est toujours plus rapide même lorsque vous augmentez sa précision avec une itération Newton-Raphson.

Intel et AMD recommandent tous deux cette technique dans leurs manuels d'optimisation. Dans les applications qui ne nécessitent pas la conformité IEEE-754, la seule raison d'utiliser div / sqrt est la lisibilité du code.

Prise de bec
la source
1
Broadwell et les versions ultérieures ont de meilleures performances de division FP, donc des compilateurs comme clang choisissent de ne pas utiliser réciproque + Newton pour scalaire sur les processeurs récents, car ce n'est généralement pas plus rapide. Dans la plupart des boucles, ce divn'est pas la seule opération, donc le débit total uop est souvent le goulot d'étranglement même en cas de divpsou divss. Voir Division en virgule flottante vs multiplication en virgule flottante , où ma réponse contient une section expliquant pourquoi le rcppsdébit n'est plus gagnant. (Ou une victoire de latence), et les chiffres divisent le débit / la latence.
Peter Cordes
Si vos exigences de précision sont si faibles que vous pouvez sauter une itération de Newton, alors oui a * rcpss(b)peut être plus rapide, mais c'est toujours plus que a/b!
Peter Cordes
5

Au lieu de fournir une réponse, cela pourrait en fait être incorrect (je ne vais pas non plus vérifier ni discuter sur le cache et d'autres choses, disons qu'ils sont identiques), je vais essayer de vous indiquer la source qui peut répondre à votre question.
La différence pourrait résider dans la façon dont sqrt et rsqrt sont calculés. Vous pouvez en savoir plus ici http://www.intel.com/products/processor/manuals/ . Je suggérerais de commencer par lire sur les fonctions du processeur que vous utilisez, il y a quelques informations, en particulier sur rsqrt (le processeur utilise une table de recherche interne avec une énorme approximation, ce qui simplifie grandement l'obtention du résultat). Il peut sembler que rsqrt est tellement plus rapide que sqrt, qu'une opération multiple supplémentaire (qui n'est pas trop coûteuse) pourrait ne pas changer la situation ici.

Edit: Quelques faits qui méritent d'être mentionnés:
1. Une fois, j'ai fait quelques micro optimalisations pour ma bibliothèque graphique et j'ai utilisé rsqrt pour calculer la longueur des vecteurs. (au lieu de sqrt, j'ai multiplié ma somme de carré par rsqrt, ce qui est exactement ce que vous avez fait dans vos tests), et cela a mieux fonctionné.
2. Calculer rsqrt en utilisant une simple table de recherche pourrait être plus facile, comme pour rsqrt, lorsque x va à l'infini, 1 / sqrt (x) va à 0, donc pour les petits x les valeurs de la fonction ne changent pas (beaucoup), alors que pour sqrt - il va à l'infini, c'est donc ce cas simple;).

Aussi, clarification: je ne sais pas où je l'ai trouvé dans les livres que j'ai liés, mais je suis à peu près sûr d'avoir lu que rsqrt utilise une table de recherche, et elle ne devrait être utilisée que lorsque le résultat n'a pas besoin d'être exact, bien que - je me trompe peut-être aussi, car c'était il y a quelque temps :).

Marcin Deptuła
la source
4

Newton-Raphson converge vers le zéro de l' f(x)utilisation des incréments égale à -f/f'f' est la dérivée.

Pour x=sqrt(y), vous pouvez essayer de résoudre f(x) = 0pour xutiliserf(x) = x^2 - y ;

Alors l'incrément est: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x qui a une division lente en elle.

Vous pouvez essayer d'autres fonctions (comme f(x) = 1/y - 1/x^2) mais elles seront tout aussi compliquées.

Regardons 1/sqrt(y)maintenant. Tu peux essayerf(x) = x^2 - 1/y , mais ce sera tout aussi compliqué: dx = 2xy / (y*x^2 - 1)par exemple. Un autre choix non évident pour f(x)est:f(x) = y - 1/x^2

Ensuite: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

Ah! Ce n'est pas une expression triviale, mais vous n'avez que des multiplications, pas de division. => Plus rapide!

Et: l'étape de mise à jour complète new_x = x + dx se lit alors:

x *= 3/2 - y/2 * x * x ce qui est facile aussi.

skal
la source
2

Il y a déjà un certain nombre d'autres réponses à cela d'il y a quelques années. Voici ce que le consensus a réussi:

  • Les instructions rsqrt * calculent une approximation de la racine carrée réciproque, bonne à environ 11-12 bits.
  • Il est implémenté avec une table de recherche (c'est-à-dire une ROM) indexée par la mantisse. (En fait, il s'agit d'une table de consultation compressée, similaire aux tables mathématiques d'autrefois, utilisant des ajustements aux bits de poids faible pour économiser sur les transistors.)
  • La raison pour laquelle il est disponible est qu'il s'agit de l'estimation initiale utilisée par le FPU pour l'algorithme de racine carrée «réelle».
  • Il existe également une instruction réciproque approximative, rcp. Ces deux instructions sont un indice sur la façon dont le FPU implémente la racine carrée et la division.

Voici ce que le consensus a mal tourné:

  • Les FPU de l'ère SSE n'utilisent pas Newton-Raphson pour calculer les racines carrées. C'est une excellente méthode dans le logiciel, mais ce serait une erreur de l'implémenter de cette façon dans le matériel.

L'algorithme NR pour calculer la racine carrée réciproque a cette étape de mise à jour, comme d'autres l'ont noté:

x' = 0.5 * x * (3 - n*x*x);

C'est beaucoup de multiplications dépendant des données et une soustraction.

Ce qui suit est l'algorithme que les FPU modernes utilisent réellement.

Étant donné b[0] = n, supposons que nous puissions trouver une série de nombres Y[i]tels que l' b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2approche 1. Considérons ensuite:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

x[n]Approches sqrt(n)et y[n]approches clairement 1/sqrt(n).

Nous pouvons utiliser l'étape de mise à jour Newton-Raphson pour la racine carrée réciproque pour obtenir un bon Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Ensuite:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

et:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

La prochaine observation clé est la suivante b[i] = x[i-1] * y[i-1]. Alors:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Ensuite:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

Autrement dit, étant donné les x et y initiaux, nous pouvons utiliser l'étape de mise à jour suivante:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

Ou, même plus chic, nous pouvons définir h = 0.5 * y. C'est l'initialisation:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

Et voici l'étape de mise à jour:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

C'est l'algorithme de Goldschmidt, et il a un énorme avantage si vous l'implémentez dans le matériel: la "boucle interne" est constituée de trois multi-ajouts et rien d'autre, et deux d'entre eux sont indépendants et peuvent être pipelined.

En 1999, les FPU avaient déjà besoin d'un circuit d'ajout / soustraction en pipeline et d'un circuit de multiplication en pipeline, sinon le SSE ne serait pas très "streaming". Un seul de chaque circuit était nécessaire en 1999 pour implémenter cette boucle interne de manière entièrement pipeline sans gaspiller beaucoup de matériel uniquement sur la racine carrée.

Aujourd'hui, bien sûr, nous avons fusionné multiply-add exposé au programmeur. Encore une fois, la boucle interne est constituée de trois FMA en pipeline, qui sont (encore) généralement utiles même si vous ne calculez pas les racines carrées.

Pseudonyme
la source
1
En relation: Comment fonctionne sqrt () de GCC après sa compilation? Quelle méthode de racine est utilisée? Newton-Raphson? a des liens vers les conceptions d'unités d'exécution matérielles div / sqrt. Rsqrt vectorisé rapide et réciproque avec SSE / AVX en fonction de la précision - une itération de Newton dans le logiciel, avec ou sans FMA, à utiliser avec _mm256_rsqrt_ps, avec l'analyse des performances Haswell. Ce n'est généralement une bonne idée que si vous n'avez pas d'autre travail dans la boucle et que vous gêneriez fortement le débit du diviseur. HW sqrt est unique et est donc bien mélangé avec d'autres travaux.
Peter Cordes
-2

C'est plus rapide parce que ces instructions ignorent les modes d'arrondi, et ne gèrent pas les exceptions de virgule flottante ou les nombres dernormalisés. Pour ces raisons, il est beaucoup plus facile de canaliser, de spéculer et d'exécuter d'autres instructions fp dans le désordre.

Witek
la source
Manifestement faux. FMA dépend du mode d'arrondi actuel, mais a un débit de deux par horloge sur Haswell et les versions ultérieures. Avec deux unités FMA entièrement pipelinées, Haswell peut avoir jusqu'à 10 FMA en vol à la fois. La bonne réponse est rsqrtla précision beaucoup plus faible, ce qui signifie beaucoup moins de travail à faire (ou pas du tout?) Après une recherche de table pour obtenir une estimation de départ.
Peter Cordes