Extraction de racine carrée

La méthode de Héron est une méthode historique développée par les Babyloniens. En termes plus modernes, c'est un cas particulier de la méthode de Newton.

Pour déterminer la racine carrée du nombre (positif) a, il convient dès lors de considérer la suite définie par récurrence de la façon suivante : ${\displaystyle \forall n\in \mathbb {N} \quad x_{n+1}={\frac {x_{n}+{\frac {a}{x_{n}}}}{2}},}$ de premier terme x₀ > 0 choisi si possible « assez proche » de √a, en général la partie entière de √a.

La vitesse de convergence des approximations successives vers la valeur exacte est quadratique.

On trouve dans un manuscrit indien, dit manuscrit de Bakshali, datant peut-être du VII^e siècle, une correction différente de la méthode de Héron, la nouvelle valeur approchée étant ${\displaystyle r+s-{\frac {s^{2}}{2(r+s)}}}$. Elle est équivalente à appliquer deux fois de suite la méthode de Héron[1]. L'itération de cette dernière méthode donne une vitesse de convergence bi-quadratique :

${\displaystyle {\sqrt {a^{2}+\gamma }}=a+{\frac {\gamma }{2a}}-{\frac {(\gamma /2a)^{2}}{2(a+\gamma /2a)}}.}$

Considérons les suites u et v définies par :

• ${\displaystyle u_{0}=1,v_{0}=a}$,
• ${\displaystyle u_{n+1}={\frac {2}{{\frac {1}{u_{n}}}+{\frac {1}{v_{n}}}}}}$ est la moyenne harmonique de u_n et v_n,
• ${\displaystyle v_{n+1}={\frac {u_{n}+v_{n}}{2}}}$ est la moyenne arithmétique de u_n et v_n.

Les suites u et v sont adjacentes, et convergent vers la même limite : √a. L'erreur est majorée par la différence v_n - u_n. La suite v n'est autre que la suite obtenue en itérant la méthode de Héron à partir de la valeur 1.

Remarquons l'originalité de cette présentation qui mêle moyennes harmonique, géométrique et arithmétique. En effet, √a n'est autre que la moyenne géométrique de 1 et de a, et si l'on remplace ${\displaystyle u_{0}}$ par un réel strictement positif quelconque b, les suites u et v convergent vers la moyenne géométrique √ab de a et b.

L'introduction de la notation décimale des nombres par position a permis de développer un algorithme tirant parti de cette notation. On en trouve mention dans un ouvrage du mathématicien indien Âryabhata, vers 499 apr. J.-C.[2] Il a été utilisé pendant plusieurs siècles et jusqu'au milieu du XX^e siècle, avant l'invention des calculateurs électroniques. Âryabhata donne également un algorithme comparable permettant de calculer des racines cubiques.

On sépare les chiffres du nombre par paires en commençant à partir de la virgule. On place le nombre dont on veut extraire la racine en haut, de la même façon que lorsqu'on effectue une division selon la méthode classique ; la racine carrée sera inscrite à la place réservée normalement au diviseur dans une division posée classique.

À chaque étape :

• on abaisse la paire de chiffres la plus significative non encore utilisée et on la place au côté d’un reste éventuel de l'étape précédente (initialement nul) ;
• soit r le résultat intermédiaire de la racine carrée obtenu précédemment (égal à zéro au début). On cherche le plus grand chiffre x tel que le nombre y = (20r + x)x ne dépasse pas le reste courant ;
• on complète r en plaçant la décimale x à sa droite, pour former le nouveau résultat intermédiaire ;
• on soustrait y de la valeur courante pour former le nouveau reste ;
• si le reste est nul et qu’il n’y a plus de chiffre à abaisser alors l’algorithme se termine sinon on recommence.

On remarque que la recherche de x en négligeant le terme en x² n'est autre que la méthode de Héron.

Jusqu’au XX^e siècle on utilisait couramment cet algorithme en accélérant les calculs à l’aide d’un abaque formé d’un jeu de réglettes : les bâtons de Napier.

Bien que décrite ici pour des nombres écrits en base 10, la procédure fonctionne dans n’importe quelle base de numération, par exemple la base 2.

On pose le problème comme la détermination de la racine positive du polynôme P(X) = X² - a. Soit n le plus grand entier tel que P(n) soit négatif ou nul. La racine de a est un réel compris entre n et n + 1. On pose alors X = n + Y/10, dont on déduit P(X) = P(n + Y/10) = P₁(Y). La racine carrée de a est alors égale à n + r/10 où r est racine du polynôme P₁.

Si m est le plus grand entier tel que P₁(m) est négatif, la racine de a est alors comprise entre n + m/10 et n + b/10 + 1/10. On pose alors Y = b + Z/10, d'où P₁(Y) = P₂(Z), puis on procède par récurrence.

La méthode de Ruffini-Horner permet d'effectuer facilement les changements de variables.

Cette méthode enseignée parfois au XIX^e siècle a l'avantage de se résumer à une série d'additions (au maximum 9 pour chaque décimale cherchée) d'impairs consécutifs[3].

Elle consiste à exploiter la table des carrés successifs et de leur différence, en remarquant que, pour tout entier n, (n+1)² - n² = 2n+1. Balayer la table des carrés consiste donc à ajouter des impairs successifs. Après avoir découpé le nombre en tranches de deux chiffres en commençant par la droite, on recherche le carré immédiatement inférieur au groupe le plus à gauche. Soit n, la racine carrée de ce carré immédiatement inférieur. On multiple l'entier n trouvé par 10 et on balaie la table des carrés à partir de 10n (et 100n²) pour s'approcher du nombre formé des deux groupes les plus à gauche. Ce nombre n₁ étant trouvé, on le multiplie par 10 pour parcourir la table des carrés à partir de 10n₁, etc.

C'est ce même principe qui est utilisé dans l'extraction de racine carrée par la méthode du goutte à goutte.

On peut raccourcir l'exécution de l'algorithme : avant de balayer la table des carrés à partir de 10n, la facile calculabilité des carrés d'entiers se terminant par 5 peut permettre, dans certains cas, de s'affranchir de cinq additions en testant si (10n+5)², soit 100n(n+1) + 25, reste aussi inférieure (ou pas) au nombre dont on cherche la racine. Si c'est le cas, il vaut mieux poursuivre l'algorithme à partir de 10n+5 que de 10n.

Par exemple, on cherche l'arrondi au dixième près de la racine carrée de 4700. Pour respecter cette précision, on est amené à chercher la racine de 4700 × 100 = 470000. Il suffira alors de diviser la racine finale par √100 = 10.

Comme 6² = 36 < 47 < 49 = 7², on a n = 6. Mais 47 étant plus proche de 49 que de 36, √47 sera plus proche de 7 que de 6.

Au lieu de commencer à balayer à 10n = 60 et 100 n² = 3600, on peut, dans ce cas, initier le processus à 10n+5 = 65 : 65² = 6 × 7 × 100 + 25 = 4225. On a bien, comme prévu, 4225 < 4700.

L'algorithme (détaillé dans l'exemple ci-dessus) amène, petit à petit, à n2
4 = 685² = 469225 auquel on ajoute 2 n₄ + 1 = 1370, d'où n2
5 = 686² = 470595 qui est plus proche de 470000 que n2
4. On en déduit : ${\displaystyle {\sqrt {4700}}\approx 68,6}$ au dixième près.

Cette méthode[4] très simple a la particularité de n'utiliser que des opérations très simples : addition, soustraction et ajout de 0. Considérons les suites a et b définies par :

• ${\displaystyle a_{0}=5m,b_{0}=5}$,
• Si ${\displaystyle a_{n}\geq b_{n}}$ alors ${\displaystyle a_{n+1}=a_{n}-b_{n}}$ et ${\displaystyle b_{n+1}=b_{n}+10}$
• Sinon, ${\displaystyle a_{n+1}=100a_{n}}$ (ce qui revient à rajouter deux zéros à la fin de ${\displaystyle a_{n}}$) et ${\displaystyle b_{n+1}=10(b_{n}-5)+5}$ (ce qui revient à insérer un zéro juste avant le dernier chiffre de ${\displaystyle b_{n}}$ car cette dernière termine toujours par un 5)

Ainsi, les chiffres de ${\displaystyle b_{n}}$ approchent de plus en plus les chiffres de ${\displaystyle {\sqrt {m}}}$. À noter que ${\displaystyle b_{n}}$ reste un entier (de plus en plus grand) donc ce n'est pas ${\displaystyle b_{n}}$ qui tend vers ${\displaystyle {\sqrt {m}}}$ mais seulement ses chiffres dans sa représentation décimale.

Exemple :√3 ≈ 1,732 par additions et soustractions

${\displaystyle b_{0}=5}$, ${\displaystyle a_{0}=15}$

${\displaystyle n}$	${\displaystyle a_{n}}$	${\displaystyle b_{n}}$	Commentaires
0	15	5	On peut soustraire
1	10	15	On ne peut plus soustraire
2	1000	105	On peut soustraire
3	895	115	On peut soustraire
4	780	125	On peut soustraire
5	655	135	On peut soustraire
6	520	145	On peut soustraire
7	375	155	On peut soustraire
8	220	165	On peut soustraire
9	55	175	On ne peut plus soustraire
10	5500	1705	On peut soustraire
11	3795	1715	On peut soustraire
12	2080	1725	On peut soustraire
13	355	1735	On ne peut plus soustraire
14	35500	17305	On peut soustraire
15	18195	17315	On peut soustraire
16	880	17325	on ne peut plus soustraire
...

Si, au lieu de prendre 5m, on en prend des troncatures par tranches de deux chiffres et, au lieu d'ajouter à a_n deux zéros, on abaisse les tranches suivantes, on aboutit à la variante par 5 de la méthode par goutte à goutte.

La fraction continue d'un irrationnel est la suite de ses approximations « optimales » par des fractions, c'est-à-dire que si p/q est une des valeurs de cette suite, alors aucune fraction de a/b avec b < q n'approche plus précisément le réel. Dans le cas particulier des racines carrées, on calcule relativement simplement cette suite, ainsi qu'une sous-suite qui correspond à un cas particulier de la méthode de Héron.

On peut également procéder par dichotomie à condition de disposer d'un encadrement de la racine carrée cherchée. On peut pour cela utiliser le nombre de chiffres du nombre dont on cherche la racine carrée. Cette racine aurait moitié moins de chiffres. Ainsi, si le nombre possède 1 024 chiffres, sa racine carrée en possèdera entre 511 et 513. On peut également utiliser d'abord les méthodes précédentes pour obtenir une première valeur approchée de la racine carrée avant de procéder à la dichotomie.

L'algorithme de dichotomie est le suivant. Il évite de procéder à des divisions (autre que la division par 2 qui n'est qu'un décalage de registre si les nombres sont codés en binaire. Cette division est notée shr 1 ci-dessous).

function Racine_64(C: int64): int64;
var
 A, B, D, D2: int64;
begin
  A := borne basse;
  B := borne haute;
  repeat
    D := (A + B) shr 1;
    D2 := D * D; ⇐ on élève au carré avant de tester
    if D2 > C then        
      A := D - 1
    else
      if C > D2 then
        B := D + 1
      else
        Result := A;
  until B > A;
end;

La même méthode s'applique pour les racines n-ièmes, en remplaçant le calcul de D2 = D*D par le calcul de D^n.

La méthode de dichotomie a cependant une vitesse de convergence plus lente que l'itération de la méthode de Héron.

La racine carrée d'un nombre complexe Z = a + ib sera un nombre complexe z = x + iy tel que :

${\displaystyle a=x^{2}-y^{2}\,b=2xy}$

En posant m² = |z|² = x² + y², le carré du module de Z, on peut établir une formule générale :

${\displaystyle z={\begin{cases}\pm {\frac {1}{2}}\left[{\sqrt {2(m+a)}}+\mathrm {i} {\sqrt {2(m-a)}}\right]&{\textrm {si}}\ b>0\\\pm {\frac {1}{2}}\left[{\sqrt {2(m+a)}}-\mathrm {i} {\sqrt {2(m-a)}}\right]&{\textrm {si}}\ b<0\end{cases}}}$

Le calcul revient donc à extraire trois racines carrées : le calcul de m, puis les racines carrées de m + a et m - a.

On cherche à résoudre l'équation x² ≡ a mod n. On distingue alors trois cas[5]:

• n est un nombre premier (par exemple 19) ;
• n est une puissance de nombre premier (par exemple 81 = 3⁴) ;
• n est un nombre composé (par exemple 35 = 5×7).

S'il existe une solution, c'est-à-dire si a est un résidu quadratique modulo n, on dispose d'algorithmes efficaces pour la trouver dans les deux premiers cas. Si n est un nombre composé, on ne dispose à l'heure actuelle d'un algorithme efficace que si la factorisation de n est connue[5]. De plus on peut montrer que s'il existait un algorithme efficace pour calculer une racine carrée sans disposer de la factorisation de n, cet algorithme pourrait être utilisé pour obtenir un algorithme efficace de factorisation. On ne pourra donc pas calculer efficacement de racines carrées modulo n tant que l'on ne pourra pas factoriser efficacement n'importe quel entier et réciproquement[5].