racine carrée

Bonjour,

J’essaye de reconstituer la fonction sqrt du module math par mes propres moyens pour comprendre comment elle fonctionne.

Voici mon modèle :

def sqrt(x):
    cf = 20

    resultat = 0

    e  = nb_chiffres(x)

    while True: #simule une boucle do...while
        e -= 1

        print(e)
        i = 0

        while (resultat + (i+1)*10**e)**2 < x and i <= 9:
            i += 1


        resultat += i*10**e

        if resultat + i*10**e == resultat and i != 0 or e == -20: #soit on a saturé l'espace mémoire alloué au nombre, soit on atteint une sécurité
            break
    return resultat

Il fonctionne plutôt bien, mais il y a quelque chose qui me chiffonne…

Quand j’exécute ma fonction avec 2 en argument, j’obtiens 1.414213562373095, mais la fonction du module math retourne 1.4142135623730951. En creusant un peu, j’ai découvert qu’en fait, le nombre trouvé par mon algorithme était 1.414213562373094999…, mais que Python arrondissait parce qu’il n’avait plus la place.

Je n’arrive pas à comprendre pourquoi il réagit comme ça…

Merci d’avance pour votre aide,

@flopy78

En effet les nombres flottants disposent d’un espace mémoire limité et donc nécessitent des arrondis, voir cet article à ce propos. Donc ce n’est pas Python qui arrondit, c’est la norme des nombres flottants qui veut ça. math.sqrt(2) renvoie la meilleure approximation possible en tenant compte de ça.

Comme ton code semble calculer le résultat décimale par décimale, tu ne retombes pas sur le même arrondi que celui du nombre renvoyé par sqrt.

Est-ce grave ? Pas vraiment : dans tous les cas tu n’obtiendras qu’un résultat approché puisque la racine de 2 est irrationnelle. Il diffère du résultat de sqrt mais semble assez valide (juste un poil moins proche).

>>> from decimal import Decimal
>>> abs(2 - Decimal(math.sqrt(2))**2)
Decimal('2.73432346306E-16')
>>> abs(2 - Decimal(sqrt(2))**2)
Decimal('3.54604637167E-16')

Les deux nombres ne sont pas égaux mais très proches. Et étant donné qu’il est déconseillé de tester l’égalité stricte entre deux flottants (à part si c’est précisément ce que l’on veut faire mais ce n’est généralement pas le cas) en raison des problèmes d’arrondis évoqués plus haut qui mènent à des erreurs de calcul, la proximité est tout ce qui devrait t’intéresser.

La fonction isclose du module math te permet de tester cela, et de préciser la tolérance.

>>> import math
>>> math.isclose(math.sqrt(2), sqrt(2))
True

En revanche ton algorithme m’a l’air plutôt inefficace : comme tu calcules décimale par décimale c’est assez lent à converger (et la boucle intérieure ralentit encore les choses). Une idée serait de t’orienter vers la méthode de Héron/Newton.

Je déplace vers le bon forum.

La méthode de Héron-Newton suppose qu’on part d’une première approximation.

Le petit code suivant illustre ce qui se passe si on choisit mal cette approximation.

def sqrt2(x, a):
    c = 0
    x0 = a
    while True:
        x1 = (x/x0 + x0) / 2
        if abs(x0-x1) < 1e-15: break
        c += 1
        x0 = x1
    print(c)
    return x0
print(sqrt2(1000000, 2))
print(sqrt2(1000000, 707))
print(sqrt2(1000000, 1414))

qui donne:

14                                                                                                                      
1000.0                                                                                                                  
5                                                                                                                       
1000.0                                                                                                                  
5                                                                                                                       
1000.0                                                                                                                   

Dans les ordinateurs, il existe des instructions pour extraire la mantisse et l’exposant des nombres en point flottant. Une approximation facile à réaliser en assembler est d’associer la moitié de l’exposant avec la mantisse initiale. C’est moins évident à faire en Python.

Ce qui reste tout à fait raisonnable, une dizaine d’itérations pour converger en partant d’une valeur éloignée du résultat.

Si l’on compare avec le code plus haut, pour 1000000 il converge en 21 itérations de la boucle principale (149 si l’on compte la boucle intérieure).

Salut,

La plupart des opérations sur les nombres flottants ne sont pas implémentées par les gens de chez Python, ils se contentent d’appeler les fonctions offertes par la libc. Il y a de très grandes chances pour que math.sqrt appelle ultimement une instruction processeur du genre sqrtsd. Celle-ci est implémentée, au niveau hardware, avec une méthode similaire à ce que tu as implémentée, mais directement en base 2 et en calculant plusieurs chiffres à la fois (voir ce papier par exemple). Les autres intervenants parlent beaucoup performance, mais à mon avis passent à côté de ce qui est important dans cette discussion. Si ce qui t’intéresse est savoir comment sqrt arrive à être rapide, c’est un problème d’architecture processeur. Il faut aller voir les implémentations d’instructions comme sqrtsd et sqrtpd, et donc ce qui marche bien pour le matériel (et comme indiqué dans le lien plus haut, c’est pas forcément ce qui marche bien quand on manipule les flottants à plus haut niveau dans du software). Si ce qui t’intéresse est savoir comment on peut développer des méthodes diverses pour calculer des racines carrées, ça devient un problème d’analyse numérique. Je t’invite à consulter cette page.

C’est un peu étrange de dire ça alors que mon message parle essentiellement de la représentation des nombres flottants et de ce qui fait que les deux résultats diffèrent bien que son calcul soit juste.
Ce n’est que dans les deux dernières phrases que j’évoque accessoirement la vitesse de convergence.