La Racine Carrée Inverse Extrêmement Rapide et Insolite de John Carmack (Quake III)

Question

John Carmack a une fonction spéciale dans le code source de Quake III qui calcule la racine carrée inverse d'un nombre à virgule flottante, 4 fois plus rapidement que la version régulière (float)(1.0/sqrt(x)), incluant une constante étrange 0x5f3759df. Voir le code ci-dessous. Quelqu'un peut-il expliquer ligne par ligne ce qui se passe exactement ici et pourquoi cela fonctionne beaucoup plus rapidement que l'implémentation régulière ?

float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;
  i  = 0x5f3759df - ( i >> 1 );
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) );

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) );
  #endif
  #endif
  return y;
}

Answer 1

À titre d'information. Carmack ne l'a pas écrit. Terje Mathisen et Gary Tarolli se sont tous deux partiellement (et très modestement) crédités pour cela, ainsi que pour certains autres crédits.

Comment la constante mythique a été dérivée reste quelque peu un mystère.

Pour citer Gary Tarolli:

Ce qui effectivement exécute un calcul en virgule flottante en entier - il a fallu longtemps pour comprendre comment et pourquoi ceci fonctionne, et je ne me souviens plus des détails.

Une constante légèrement meilleure, développée par un mathématicien expert (Chris Lomont) essayant de comprendre comment fonctionnait l'algorithme original est :

float InvSqrt(float x)
{
    float xhalf = 0,5f * x;
    int i = *(int*)&x;              // obtenir les bits pour la valeur flottante
    i = 0x5f375a86 - (i >> 1);      // donne un essai initial y0
    x = *(float*)&i;                // convertir les bits en float
    x = x * (1,5f - xhalf * x * x); // étape de Newton, répéter augmente la précision
    return x;
}

Malgré cela, sa tentative initiale d'une version mathématiquement 'supérieure' de la racine carrée de id (qui a abouti à une constante presque identique) s'est révélée inférieure à celle initialement développée par Gary malgré le fait d'être mathématiquement beaucoup plus 'pure'. Il ne pouvait pas expliquer pourquoi celle de id était si excellente, de mémoire.

Answer 2

Bien sûr de nos jours, il s'avère beaucoup plus lent que d'utiliser simplement la racine carrée d'un FPU (surtout sur 360/PS3), car le passage entre les registres de nombres flottants et entiers induit un load-hit-store, tandis que l'unité de calcul en virgule flottante peut faire la racine carrée réciproque matériellement.

Cela montre simplement comment les optimisations doivent évoluer au fur et à mesure que la nature du matériel sous-jacent change.

Answer 3

Greg Hewgill et IllidanS4 ont donné un lien avec une excellente explication mathématique. Je vais essayer de résumer ici pour ceux qui ne veulent pas entrer trop dans les détails.

Toute fonction mathématique, à quelques exceptions près, peut être représentée par une somme de polynômes :

y = f(x)

peut être exactement transformé en :

y = a0 + a1*x + a2*(x^2) + a3*(x^3) + a4*(x^4) + ...

Où a0, a1, a2,... sont des constants. Le problème est que pour de nombreuses fonctions, comme la racine carrée, pour une valeur exacte cette somme a un nombre infini de membres, elle ne se termine pas à un certain x^n. Mais, si nous nous arrêtons à un certain x^n, nous aurions quand même un résultat jusqu'à une certaine précision.

Donc, si nous avons :

y = 1/sqrt(x)

Dans ce cas particulier, ils ont décidé de se débarrasser de tous les membres polynomiaux au-dessus de deux, probablement à cause de la vitesse de calcul :

y = a0 + a1*x + [...discarded...]

Et la tâche est maintenant de calculer a0 et a1 afin que y ait la moins de différence possible par rapport à la valeur exacte. Ils ont calculé que les valeurs les plus appropriées sont :

a0 = 0x5f375a86
a1 = -0.5

Donc quand vous mettez cela dans l'équation vous obtenez :

y = 0x5f375a86 - 0.5*x

Ce qui est la même que la ligne que vous voyez dans le code :

i = 0x5f375a86 - (i >> 1);

Edit : en fait ici y = 0x5f375a86 - 0.5*x n'est pas la même que i = 0x5f375a86 - (i >> 1); car le décalage du flottant en tant qu'entier divise non seulement par deux mais divise aussi l'exposant par deux et provoque d'autres artefacts, mais cela revient toujours à calculer certains coefficients a0, a1, a2... .

À ce stade, ils ont découvert que la précision de ce résultat n'était pas suffisante pour l'objectif. Ils ont donc effectué une seule étape de l'itération de Newton pour améliorer la précision du résultat :

x = x * (1.5f - xhalf * x * x)

ils auraient pu faire d'autres itérations dans une boucle, chacune améliorant le résultat, jusqu'à ce que la précision requise soit atteinte. C'est exactement ainsi que cela fonctionne dans le CPU/FPU ! Mais il semblerait qu'une seule itération ait été suffisante, ce qui était également une bénédiction pour la vitesse. Le CPU/FPU effectue autant d'itérations que nécessaire pour atteindre la précision du nombre à virgule flottante dans lequel le résultat est stocké et il dispose d'un algorithme plus général qui fonctionne pour tous les cas.

---

En bref, ce qu'ils ont fait c'est :

Utiliser (presque) le même algorithme que le CPU/FPU, exploiter l'amélioration des conditions initiales pour le cas spécial de 1/sqrt(x) et ne pas calculer jusqu'à la précision à laquelle le CPU/FPU ira mais s'arrêter avant, ce qui permet de gagner en vitesse de calcul.

Answer 4

J'étais curieux de voir quelle était la constante en tant que float, donc j'ai simplement écrit ce bout de code et j'ai googlé l'entier qui en est sorti.

long i = 0x5F3759DF;
float* fp = (float*)&i;
printf("(2^127)^(1/2) = %f\n", *fp);
//Output
//(2^127)^(1/2) = 13211836172961054720.000000

Il semble que la constante est "une approximation entière de la racine carrée de 2^127 mieux connue sous la forme hexadécimale de sa représentation en virgule flottante, 0x5f3759df" https://mrob.com/pub/math/numbers-18.html

Sur le même site, tout est expliqué. https://mrob.com/pub/math/numbers-16.html#le009_16

Answer 5

Conformément à cet article intéressant écrit il y a quelque temps...

La magie du code, même si vous ne pouvez pas le suivre, se manifeste avec la ligne i = 0x5f3759df - (i>>1);. Simplifié, Newton-Raphson est une approximation qui commence par une supposition et la affine avec une itération. Profitant de la nature des processeurs x86 32 bits, i, un entier, est initialement défini à la valeur du nombre en virgule flottante dont vous voulez inverser le carré, en utilisant une conversion en entier. i est ensuite défini à 0x5f3759df, moins lui-même décalé d'un bit vers la droite. Le décalage à droite supprime le bit de poids le moins élevé de i, le divisant essentiellement par deux.

C'est vraiment intéressant à lire. Ceci n'est qu'une petite partie.