La plupart des processeurs modernes fournissent une opération de racine carrée de manière native. Donc, avec un langage de programmation typique sur du matériel moderne typique, c'est presque certainement ce que l'opération sera en fin de compte.
Parlons donc de la méthode que presque tous les CPU modernes polyvalents utilisent.
Représentation des nombres
La première chose que nous devons comprendre est la façon dont un ordinateur représente réellement le genre de nombres dont nous parlons. Ils sont stockés dans une sorte de notation scientifique.
Maintenant, le type de notation scientifique que vous connaissez peut-être où un nombre comme -134,5 est représenté comme [math]-1,345 times 10^{2}[/math].
Tout nombre qui n'est pas zéro peut être représenté sous cette forme, qui comporte conceptuellement trois parties : le signe (c'est-à-dire que le nombre est positif ou négatif), la mantisse (un nombre compris entre 1 et 10, incluant 1 mais pas 10, dans ce cas 1.345), et l'exposant (la puissance à laquelle le radix est élevé, dans ce cas 2 ; notez que l'exposant peut être négatif ou nul en général).
La représentation que presque tous les ordinateurs modernes utilisent est très similaire, sauf que le nombre est stocké en utilisant une représentation binaire plutôt que décimale. Ainsi, en ignorant le bit de signe et le zéro (puisque trouver la racine carrée de zéro est trivial, et que la racine carrée d'un nombre négatif n'est pas un nombre réel), les nombres positifs sont en fait stockés sous la forme [math]m times 2^{e}[/math] où m est un nombre compris entre 1 et 2 (y compris 1 mais pas 2).
Cette représentation est appelée "virgule flottante", en raison de la façon dont vous pouvez déplacer le point binaire (analogue au point décimal que vous connaissez peut-être) en ajustant l'exposant.
La raison pour laquelle tout cela est important est que le CPU utilise cette représentation pour aider à calculer les racines carrées. Supposons que l'exposant soit pair. Alors :
[math]sqrt{m times 2^{2epsilon}} = sqrt{m} times 2^{epsilon}[/math]
De même, si l'exposant est impair, alors:
[math]sqrt{m times 2^{2epsilon+1} = sqrt{2} sqrt{m} times 2^{epsilon}[/math]
Rappelez-vous que m est dans l'intervalle [1,2). Cela signifie que nous avons réduit le problème du calcul de la racine carrée de n'importe quel nombre à celui du calcul de la racine carrée d'un nombre dans cet intervalle. Ou, si nous ne voulons pas précalculer [math]sqrt{2}[/math], un nombre compris dans l'intervalle [1,4]. Quoi qu'il en soit, nous avons considérablement simplifié le problème, car nous pouvons maintenant utiliser des méthodes qui se comportent bien dans cette plage de nombres.
Multiplication et division
Alors, maintenant que nous sommes arrivés jusqu'ici, vous pourriez penser (si vous avez regardé les autres réponses) que la solution consiste maintenant à utiliser une méthode comme l'itération de Newton-Raphson pour calculer la racine carrée. Pour calculer la racine carrée de n, choisissez une estimation initiale [math]x_0[/math] et itérez :
[math]x_{i+1} = frac{1}{2} left( x_i + frac{n}{x_i} right)[/math]
Cette réponse est fausse. Enfin, elle n'est pas fausse dans le sens où elle vous donnera une réponse correcte. Mais elle est fausse en ce qu'aucun concepteur de processeur qui se respecte (ou tout auteur de bibliothèque s'il devait l'implémenter en logiciel) n'implémenterait la racine carrée en virgule flottante de cette façon.
La division par deux est une opération très bon marché en binaire (qui, rappelez-vous, est en base 2, donc il s'agit juste de décaler le point binaire d'une place). Cependant, l'autre division est une opération très coûteuse, et dans ce cas, vous devez effectuer l'opération plusieurs fois.
La division est si coûteuse, en fait, que les processeurs modernes utilisent un algorithme itératif qui est similaire à (mais pas réellement) l'itération de Newton-Raphson pour effectuer la division ! Il est clair que nous ne voulons pas faire cela en matériel, dans une boucle interne.
Sur le matériel informatique moderne, il est beaucoup moins cher d'effectuer une opération de multiplication qu'une opération de division. La raison en est un peu complexe à expliquer ; elle a à voir avec le fait que nous pouvons faire tenir beaucoup plus de transistors sur une puce qu'auparavant, et la multiplication est un problème que vous pouvez déplacer sur beaucoup de transistors de manière efficace. Recherchez les arbres de Wallace si vous êtes intéressé par les détails.
En tout cas, le point est que la multiplication est une opération comparativement bon marché. Donc, si nous pouvons mettre en œuvre l'opération de racine carrée en termes de multiplication plutôt que de division, ce serait mieux.
Newton-Raphson, prise deux
Vient maintenant la première intuition clé : au lieu de calculer la racine carrée, calculez la réciproque de la racine carrée. Autrement dit, au lieu de [math]sqrt{n}[/math], calculez [math]frac{1}{sqrt{n}}[/math]. Il s'avère que c'est un nombre beaucoup plus facile à calculer, et si vous avez besoin de la racine carrée, multipliez ce nombre par n et vous avez terminé.
La méthode Newton-Raphson pour la racine carrée réciproque ressemble à ceci. Si [math]x_0[/math] est une estimation initiale à [math]frac{1}{sqrt{n}}[/math], itérer :
[math]x_{i+1} = frac{1}{2} x_i left( 3 - n {x_i}^2 right)[/math]
Encore, la division par deux est assez bon marché, et tout le reste n'est que multiplication et addition/soustraction.
C'est une excellente façon de l'implémenter dans un logiciel. En outre, il convient de souligner que dans de nombreuses situations pratiques (par exemple, la normalisation d'un vecteur en utilisant le théorème de Pythagore), la racine carrée réciproque est en fait plus utile que la racine carrée, puisque la division par une racine carrée est moins efficace que la multiplication par une racine carrée réciproque.
Cependant, ce n'est pas ainsi qu'elle est généralement implémentée dans le matériel.
Algorithme de Goldschmidt
Nous pouvons maintenant examiner l'algorithme de Goldschmidt, qui calcule la racine carrée et la racine carrée réciproque ensemble.
Donné [math]b_0 = n[/math], si nous pouvons trouver une série de nombres [math]Y_i[/math] telle que [math]b_n = b_0 {Y_0}^2 {Y_1}^2 cdots {Y_{n-1}}^2[/math] approche de 1, alors [math]y_n = {Y_0} {Y_1} cdots {Y_{n-1}}[/math] s'approchera de [math]frac{1}{sqrt{b_0}}[/math] et [math]x_n = b_0 {Y_0} {Y_1} cdots {Y_{n-1}}[/math] se rapprochera de [math]sqrt{b_0}[/math].
La série que nous utilisons est essentiellement l'étape de mise à jour de Newton-Raphson :
[math]begin{align*}b_i & = b_{i-1} {Y_{i-1}}^2 \_Y_i & = frac{1}{2}(3 - b_i)end{align*}[/math]
Et puis nous pouvons garder la trace de la racine carrée :
[math]N-gin{align*}x_0 & = n Y_0 N- x_{i} & = x_{i-1} Y_iend{align*}[/math]
Et la racine carrée réciproque :
[math]Nbgin{align*}y_0 & = Y_0 \\N y_{i} & = y_{i-1} Y_i end{align*}[/math]
Mais si nous gardons trace de [math]x_i[/math] et [math]y_i[/math], observez que [math]b_i = x_{i-1} y_{i-1}[/math]. Nous n'avons donc en fait jamais à garder trace de [math]b_i[/math]:
[math]begin{align*} Y_i & = frac{1}{2} left(3 - b_iright) \N & = 1 + frac{1}{2} left(1 - b_iright) \ & = 1 + frac{1}{2} left(1 - x_{i-1} y_{i-1}right) end{align*}[/math]
Et maintenant, nous n'avons pas non plus besoin de suivre [math]Y_i[/math]:
[math]begin{align*} x_i & = x_{i-1} left( 1 + frac{1}{2} left(1 - x_{i-1} y_{i-1}right) right) \N & = x_{i-1} + x_{i-1} frac{1}{2} left(1 - x_{i-1} y_{i-1}right) _i & = y_{i-1} left( 1 + frac{1}{2}) left(1 - x_{i-1} y_{i-1}right) right) & = y_{i-1} + y_{i-1} frac{1}{2} left(1 - x_{i-1} y_{i-1}right)end{align*}[/math]
Il y a ici des calculs redondants, que nous pouvons supprimer, en proposant l'algorithme suivant. Étant donné [math]Y[/math] comme approximation de [math]frac{1}{sqrt{n}}[/math], on définit:
[math]begin{align*} x_0 & = n Y \\end{align*}[/math]
Puis on itère :
[math]NBegin{align*} r_i & = frac{1}{2}left( 1 - x_i y_i right) N x_{i+1} & = x_i + x_i r_i N y_{i+1} & = y_i + y_i r_iNend{align*}[/math]
Même si la division par deux est bon marché, nous pouvons aussi l'éviter en gardant trace de [math]h_i = frac{1}{2} y_i[/math] au lieu de [math]y_i[/math]. C'est l'algorithme de Goldschmidt.
Supposons que [math]Y[/math] soit une approximation de [math]frac{1}{sqrt{n}}[/math]. On définit:
[math]begin{align*} x_0 & = n Y \h_0 & = frac{Y}{2}end{align*}[/math]
Puis on itère:
[math]begin{align*} r_i & = frac{1}{2} - x_i h_i \N x_{i+1} & = x_i + x_i r_i Nh_{i+1} & = h_i + h_i r_iNend{align*}[/math]
Alors que [math]x_i[/math] converge vers [math]sqrt{n}[/math] et [math]h_n[/math] converge vers [math]frac{1}{2sqrt{n}}[/math].
Mise en œuvre matérielle
Pour l'instant, tout va bien. C'est certainement un algorithme très simple. Mais pourquoi est-il meilleur ?
Les processeurs modernes disposent souvent d'un circuit rapide qui effectue une opération de multiplication-accumulation optimisée, souvent appelée multiplication-addition fusionnée, ou FMA pour faire court. Si vous avez consulté la référence aux arbres de Wallace plus tôt, vous devriez comprendre comment FMA pourrait être presque aussi efficace qu'une opération de multiplication directe.
FMA est l'une des primitives les plus utiles à avoir sous la main. Si vous devez évaluer un polynôme, par exemple :
[math]p(x) = a_0 + a_1 x + a_2 x^2 + cdots a_n x^n[/math]
vous pouvez utiliser la règle de Horner, qui est essentiellement un ensemble de FMA :
[math]N-gin{align*}p_{n-1} & = a_{n-1} + x a_n NNNNNN- p_{n-2} & = a_{n-2} + x p_{n-1} \\N & vdots N p_1 & = a_1 + x p_2 N p_0 & = a_0 + x p_1end{align*}[/math]
La boucle interne de l'algorithme de Goldschmidt est constituée de trois FMA et rien d'autre. C'est pourquoi c'est un avantage : vous n'avez besoin que d'un seul type de circuit (éventuellement un seul circuit ; notez que les deux seconds FMA sont indépendants et peuvent donc bénéficier du pipelining) plus une certaine logique de contrôle pour tout mettre en œuvre. Et c'est un circuit qui est utile dans de nombreuses autres opérations, donc vous ne gaspillez pas beaucoup de matériel uniquement pour l'opération de racine carrée.
L'avant-dernière pièce du puzzle est de savoir comment obtenir une bonne estimation initiale, et la réponse courte est que la meilleure méthode est d'utiliser une table de consultation. Même une table de taille modeste, parce que nous ne cherchons que des racines carrées dans une si petite plage, est assez efficace.
La dernière pièce du puzzle est : Comment savons-nous quand nous avons fini d'itérer ? Et la réponse à cela est que nous connaissons la précision des nombres que nous utilisons, et la qualité de l'estimation initiale. À partir de là, nous pouvons calculer à l'avance le nombre maximum d'itérations nécessaires. Donc, nous ne bouclons pas réellement en tant que tel, nous faisons juste l'itération un nombre fixe de fois.