TangenteX.com

L'importance de l'arithmétique d'un ordinateur

La perception commune de l'arithmétique d'un ordinateur

Comme beaucoup de pages de TangenteX, le besoin de rédiger cette page me vint en écoutant les discussions de mes proches ou de mes collègues ingénieurs (un physicien numéricien ne ferait jamais ce genre de remarque :-). J'entends par exemple :

J'ai aussi constaté que beaucoup d'utilisateurs de calculatrices ou d'ordinateurs n'ont aucune idée de la manière dont ces machines se représentent les nombres.

Qu'est-ce que l'arithmétique d'un ordinateur

Nos modèles mathématiques utilisés en physique manipulent des nombres réels ou complexes, assez peu souvent des entiers naturels. Lorsque nous voulons construire des algorithmes numériques et dérouler ces algorithmes sur un ordinateur, nous devons passer de l'ensemble des réels ou des complexes à un ensemble fini et discret qui soit compatible avec les capacités de nos ordinateurs à traiter les nombres, capacités qui n'ont pas l'infinitude des réels, des complexes ou des entiers.

Ces capacités sont définies par l'arithmétique de l'ordinateur, répartie subtilement entre le processeur (plus exactement son ALU - Unité Arithmétique et Logique) et le compilateur que nous utilisons. Processeurs et compilateurs ont fait beaucoup de progrès, mais il arrive encore qu'on obtienne des résultats différents pour un même algorithme de calcul en fonction du couple (processeur, compilateur).

Il importe donc que les aspirants physiciens qui utiliseront le calcul scientifique comprennent l'arithmétique de leur ordinateur et les sources éventuelles de désagréments ou pire. Pour illustrer le problème, voici un rappel de quelques catastrophes informatiques.

Certains bugs de calcul ont fait la Une de l'actualité, parfois de façon tragique. J'ai retenu ici trois cas significatifs de l'importance d'une bonne compréhension du traitement des nombres par un ordinateur. Dans ces trois cas, les erreurs étaient prévisibles et auraient pu être mises en évidence par des tests appropriés.

Erreurs de calcul et catastrophe

Le vol 501 d'Ariane

Comme toutes les fusées (les avions, les navires, les missiles aussi d'ailleurs...), Ariane 5 est dotée d'une centrale inertielle, instrument qui permet à l'appareil de se repérer dans l'espace et de suivre la trajectoire prévue. Une centrale inertielle est composée pour l'essentiel d'accéléromètres, de gyroscopes et de calculateurs, qui mesurent des accélérations et des angles sur les trois axes de l'espace et qui intègrent ces données mesurées pour calculer des vitesses et des positions. Les accéléromètres et les gyroscopes sont des capteurs qui convertissent des caractéristiques mécaniques (resp. accélération et angle) en grandeurs électriques qui sont alors numérisées, converties en nombres manipulables par un calculateur.

Et nous arrivons au coeur du problème, plus précisément sur l'information délivrée par un accéléromètre sur l'axe horizontal de la centrale. La centrale inertielle utilisée sur Ariane 5 était identique à celle utilisée sur Ariane 4. Après tout, les principes de guidage et les informations de navigation nécessaires sont les mêmes. Et la centrale inertielle d'Ariane 4 fonctionnait parfaitement !

Le premier tir d'Ariane 5 eut lieu le 4 juin 1996. 37 secondes après le décollage, la trajectoire de la fusée est sortie du nominal, la fusée a perdu un booster, ce qui déclencha son auto-destruction. Tout ça à cause d'un dépassement de buffer sur une donnée issue d'un accéléromètre horizontal. Il faut savoir que les accélérations horizontales max. d'Ariane 5 sont 5 fois plus fortes que celles d'Ariane 4 et que sa trajectoire après décollage est différente de celle d'Ariane 4, qui implique des accélérations horizontales importantes.

Sur Ariane 4, la valeur numérique max délivrée par le capteur était égale à 64, codée sur un entier court de 8 bits, largement suffisant pour coder un nombre entier égal à 64 au maximum. Le gros hic c'est que la valeur d'accélération délivrée par le capteur d'Ariane 5 lors du décollage est très supérieur, de l'ordre de 300 max. Et comme vous le savez, il n'est pas possible de coder 300 sur 8 bits (on code de 0 à 255). Il aura fallut un entier de 16 bits ou plus. Cela génère un dépassement de buffer et si l'exception n'est pas traitée, le programme se "plante". Ce qui arriva sur la centrale inertielle principale. Le calculateur de la centrale étant redondé, la centrale de secours a pris le relais et se planta immédiatement pour les mêmes raisons, le logiciel étant le même dans les deux cas.

Ce qu'il faut retenir de cet exemple :

Panne sur un croiseur Aegis de l'US Navy

Un autre cas d'invalidité des données d'entrée dans un système d'acquisition, moins pénalisant financièrement mais plus inquiétant. Le 21 septembre 1997, l'USS Yorktown CG-48, un croiseur de classe Aegis, fut victime d'une panne de son système de propulsion qui l'immobilisa, heureusement alors qu'il était à l'entrainement au large de la Virginie.

L'USS Yorktown fut le premier navire de l'US Navy à recevoir un système de gestion intégrée des systèmes de navigation, propulsion, communications et supervision des équipements, le système Smart Ship System (SSS) issu du programme SSPO (Smart Ship Program Office). Les systèmes de combats n'étaient pas concernés, mais que fait un croiseur sans propulsion ni communications ? Signalons que le but du programme SSPO était de réduire de 40 à 4 le nombre des personnels chargés de la conduite du navire...

SSS était composé de 27 PC tournant sous Windows NT 4.0, reliés entre eux et à un serveur par un LAN sur fibre optique. NT 4.0 n'était pas réputé pour sa fiabilité, sa robustesse et son endurance (il fallait rebooter fréquemment à cause des fuites mémoire du système).

Donc, lors de l'entrainement, un marin constate qu'une vanne est physiquement fermée alors que le superviseur de SSS (Smart Ship's Standard Machinery Control System) indique qu'elle est ouverte. Problème de capteur assez classique en supervision. Notre marin va donc tenter de modifier manuellement, dans la base de données d'états du système, l'état de la vanne en saisissant dans l'IHM du gestionnaire de la base de données la valeur 0 dans un attribut de reset de la vanne. Il se trouve que cet attribut était utilisé dans une division, au diviseur... Bien sur, le programme a produit une exception de division par zéro. Et le programme ne traitait pas cette exception, le programmeur ayant sans doute pensé que cet attribut ne pouvait pas prendre une valeur nulle. Le programme s'est planté et l'erreur s'est propagée à l'ensemble des programmes du SSS : panne de propulsion. Imaginez la même situation au combat...

Ce qu'il faut retenir de cet incident :

Echec de l'interception d'un SCUD par un Patriot sur Dhahran

Finissons ce rapide tour d'horizon par une erreur d'arrondi aux conséquences humaines désastreuses.

Le système de missiles anti-missiles Patriot (MIM-104F - PAC3) fut la star de la première guerre du Golfe. Sa mission était d'empêcher les SCUD (un missile balistique de fabrication soviétique) irakiens d'atteindre les troupes coalisées stationnées au Koweït et en Arabie Saoudite.

Le 25 février 1995, un SCUD fut lancé contre la base US Army de Dhahran (Arabie Saoudite). La batterie de Patriot qui assurait sa défense anti-aérienne a bien détecté le missile, l'a engagé et a tiré un missile pour l'intercepter. Ce missile a loupé sa cible de plus de 500 mètres, alors que la précision habituelle du système est d'environ un mètre (en mode ABM, la charge militaire est une charge cinétique : les missiles balistiques vont plus vite que l'onde d'explosion, ce qui rend inopérant les charges explosives déclenchées par proximité en ABM). Le SCUD provoqua la mort de 28 personnes et en blessant plus ou moins gravement 100 autres.

L'enquête montra qu'une erreur d'arithmétique du calculateur, plus précisément des arrondis cumulés, fut à l'origine de la catastrophe.

Pour comprendre ce qu'il s'est passé, il faut retenir deux données essentielles :

Dans le calculateur de tir, le temps est stocké en binaire dans un registre de 24 bits. Mais essayez de convertir 1/10 de seconde en binaire : vous allez être obligé d'arrondir le résultat. Oh pas de beaucoup : environ 9.5 10-8 seconde. Cela semble peu, mais comptez combien il y a de pas d'horloge au dixième de seconde dans 100 heures : 3.6 106. Ce qui représente une erreur d'arrondi du temps de 9.5*10-8*3.6*106 = 0.34 seconde. Le temps écoulé est sous-estimé de 0.34 s. Un SCUD en phase balistique vole à 1 676 m.s-1. En 0.34 seconde, il parcourt 570 mètres de plus que prévu par le calculateur : l'interception échoue.

Nous sommes en présence d'un cumul d'erreurs dues aux arrondis successifs. Pour l'anedocte si l'on peut dire, les israéliens, eux aussi utilisateurs de Patriot, ont prévenus les américains du problème le 11 février 1991. Le patch de correction a été transmis le 26 février 1991. La seule parade qu'avait indiqué le constructeur du Patriot, Raytheon, était de rebooter le calculateur régulièrement, mais il semble que les servants n'aient pas compris pourquoi...

Dans ce cas, il s'agit clairement un défaut de conception, aggravé par une absence de test d'endurance :

Autres erreurs du même genre

L'histoire de l'informatique est riche d'erreurs plus ou moins dramatiques du même genre, par exemple :

Que faut-il en conclure

En conclusion de ces quelques exemples :

Représentation des nombres

Ensemble des réels et ensemble des nombres sur un ordinateur

L'ensemble des réels

L'ensemble des réels possède une caractéristique très particulière : si vous choisissez un élement a de cet ensemble, il existe un élement b de ce même ensemble, tel que la distance entre a et b soit aussi petite que vous le souhaitez, et même infiniment petite (la notion de distance est définie sur l'ensemble des réels). L'ensemble des réels est continu. Autrement dit, il n'existe pas de trou entre deux réels voisins (le voisinage est aussi défini sur l'ensemble des réels).

L'ensemble des réels est infini : quelque soit un réel a, il est possible de construire un réel b, tel que b soit supérieur à a. Et quelque soit c, il est toujours possible de construire d tel que d soit inférieur à c (la relation d'ordre est défini sur l'ensemble des réels).

L'ensemble des nombres manipulables sur un ordinateur

Un ordinateur est une machine physique qui stocke et manipule des nombres. Ces nombres sont stockés dans une mémoire et manipulés par l'ALU d'un processeur.

Imaginons que je veuille effectuer l'addition de deux nombres réels a et b avec mon ordinateur. Il va falloir que je stocke ces deux nombres en mémoire et que mon processeur sache les manipuler. Un nombre réel peut être représenté sous forme décimale, sa suite de décimales étant infinie. Ce qui me pose problème ; je n'ai pas une quantité de mémoire infinie ! Et mon processeur mettrait un temps infini pour lire et traiter un nombre décimal de longueur infinie ! Il faut je trouve une représentation finie de mes nombres réels a et b qui soit compatible avec les capacités de ma mémoire et de mon processeur.

La finitude de cette représentation est une caractéristique fondamentale de l'arithmétique des nombres sur ordinateur. L'ensemble de ces nombres est un ensemble fini. Il n'a pas la puissance du continu, c'est à dire qu'il n'existe pas de bijection entre cet ensemble et l'ensemble des réels.

Ce n'est pas un ensemble continu : il est possible de construire deux nombres A et B tels que je puisse minorer leur distance, c'est à dire que je puisse exhiber une distance la plus petite entre deux nombres qui ne tend pas vers 0.

Qu'on se le dise : un ordinateur ne sait pas manipuler des nombres réels mais une représentation finie des nombres réels ! Une grande partie du travail d'un physicien numéricien est justement de tenir compte de toutes les conséquences de cette restriction.

Pour finir, une réflexion philosophique : il n'est pas certain que notre univers soit continu... Et si notre monde microscopique était discret lui aussi ? Alors tous nos modèles continues réels seraient à adapter...

Les contraintes matérielles

Au début de l'informatique basée sur des microprocesseurs, disons au début des années 80, la taille des cases mémoires dans lesquelles sont stockions nos données était de 8 bits. Ces cases mémoires pouvaient prendre deux états 0 et 1. Il est ressort que nous disposions de 28 combinaisons pour représenter un nombre. Les processeurs manipulaient aussi des mots de 8 bits (par exemple l'Intel 8080 ou le Motorola 6800 pour les anciens...).

Les progrès arrivèrent vite : la taille des mots est passée à 16 bits puis 32 bits puis 64 bits, pour les processeurs le plus courant, et il existe même des processeurs qui traitent des mots de 128 ou 256 bits.

Il n'en reste pas moins vrai que nos mots mémoires sont toujours finis, et que se pose la question de représenter un nombre réel, infini par définition, avec un nombre de combinaisons de bits fini, même si celui est grand (232 ou 264).

Représentation d'un nombre

La forme générale

Un nombre réel x peut être représenté sous la forme générale :

\( x = \displaystyle sgn(x) \beta^{e} \sum_{k \ge 1} \dfrac{d_k}{\beta^k} \)

dans laquelle :

Par exemple, le réel -0.0156 est représenté par la forme \( \displaystyle  -0.0156 = - 10^{-1} \left(\dfrac{1}{10} + \dfrac{5}{10^2} + \dfrac{6}{10^3} \right) \).

Les nombres flottants

En pratique, on utilise une forme plus compacte, et l'on introduit la notion de "virgule flottante". Les nombres à virgule flottante, que l'on appelle plus couramment chez les numériciens les nombres flottants, s'expriment sous la forme :

\( x = sgn(x) m \beta^{e} \)

dans laquelle :

Dans la base \( \beta \) choisie, la mantisse comporte p chiffres, p désigne alors la précision de la représentation du réel x.

L'exposant e varie entre \(e_{min} \le e \le e_{max} \), qui constitue le domaine des exposants. La valeur de l'exposant traduit la place de la virgule dans l'écriture décimale du réel.

La donnée de \( {\beta, p, e_{min}, e_{max} }\) constitue un système de nombres à virgule flottante, qui permet de construire une arithmétique à virgule flottante.

"machine epsilon"

Dans beaucoup d'algorithmes numériques vous trouverez une constante particulière, notée "eps" dans pratiquement tous les langages de programmation modernes. "eps", pour "epsilon", désigne la valeur de la distance de 1 au nombre flottant le plus proche dans la représentation en virgule flottante de la machine, pour un type de représentation donnée.

Pour illustrer cette notion, voyons quelle est la valeur de eps pour différents types de nombres flottants en Python. Pour un type de variable en float32, c'est à dire un nombre flottant codé sur 32 bits ( nous verrons cela après), j'utilise la commande suivante pour obtenir la valeur de eps :
print "eps en float32 : ", finfo(float32).eps
j'obtiens la valeur :
eps :  1.19209e-07

Pour des flottants codés sur 32 bits en Python, la distance entre deux nombres flottants adjacents est donc d'environ 1,2.10-7. Cette valeur traduit la discrétisation de l'ensemble des nombres flottants codés sur 32 bits.

Le type float en Python est codé sur 64 bits. Dans ce cas, la valeur de eps est de 2.22044604925e-16, soit environ 2,2.10-16. Le pas est beaucoup plus petit, ce qui est appréciable dans beaucoup d'algorithmes !

Ce qu'il faut retenir, c'est que cela n'a aucun sens de vouloir obtenir une précision de résultat plus grande que ce que nous permet la valeur de eps pour le type de variable choisi. Vouloir atteindre une précision de 10-8 sur un résultat de type float32 n'a pas de sens. Il faut changer de type et passer en float.

Enfin, Python supporte le type float128, des flottants codés sur 128 bits. Dans ce cas, la constante eps vaut 1.08420217249e-19. C'est la plus grande précision que vous pourrez atteindre en Python, sans utilisation d'algorithmes particuliers.

La représentation physique des réels sur un ordinateur

Sur les PC et stations de travail que nous utilisons quotidiennement, les nombres réels sont représentés par des nombres flottants codés sur 32 bits, on parle alors de "simple précision" ou sur 64 bits, et l'on parle alors de "double précision". La base utilisée est la base 2.

En simple précision, 24 bits dont un bit de signe servent à coder la mantisse et 8 bits sert à coder l'exposant.
En double précision, 53 bits dont un bit de signe servent à coder la mantisse et 11 bits sert à coder l'exposant.

Aujourd'hui, tout cela est normalisé, comme nous le verrons dans le chapitre suivant, mais ce ne fut pas toujours le cas ! Par exemple, sur les IBM 3090, très utilisés en calcul scientifique quand j'ai débuté ma carrière, la base était hexadécimale (16), le nombre de chiffres de la mantisse variait entre 6 et 28 et l'exposant variait entre -64 et +63. Sur les VAX, autres machines très communes dans les années 80 et 90, la base était binaire, le nombre de chiffres de la mantisse variait entre 53 et 56 et l'exposant variait entre -1023 et +2013. Sans parler du Cray 1, supercalculateur de l'époque...

Bref, il était très difficile d'écrire un code FORTRAN (par exemple) portable d'un 3090 sur un VAX (sans parler des problèmes d'OS), qui donne le même résultat d'une machine à l'autre et même d'un compilateur à l'autre.

Un autre problème fondamental se posait. Les arithmétiques variaient d'un calculateur à l'autre, et parfois ne respectaient même pas les propriétés de l'arithmétique des réels. Par exemple, sur le Cray 1, l'addition n'était pas commutative : x+y était différent de y+x ! Imaginez la joie de programmer à l'époque !

Une norme vint heureusement mettre un terme à ce bazard : la norme IEEE-754.

La norme IEEE-754

Ses objectifs

Le début des années 80 marquèrent le début de l'industrialisation du logiciel. Les problèmes de portabilité et de reproductibilité des résultats selon les machines et les compilateurs devinrent critiques.

La norme IEEE-754 ""Standard for binary floating-point arithmetic" en 1985, et sa petite soeur l'IEEE-854 "Standard for radix-independent floating-point arithmetic" en 1987, avaient pour objectif essentiel de garantir la portabilité des codes et l'indépendance de l'arithmétique des processeurs et des compilateurs. Elles poursuivaient également d'autres objectifs, parmi les principaux :

Ce que fixe la norme

La norme IEE-754 fixe principalement :

A noter que la norme IEEE-754 fut révisée et complétée en 2008.

La norme IEEE-754 fixe beaucoup de choses, mais cela ne dispense pas de réflexion lorsqu'on écrit un code de calcul. Dans certains cas, le programmeur doit déterminer le comportement de son code. Dans ce cas, un seul commandement : le code ne doit pas se "planter" ! Le traitement doit continuer.

Notons aussi que l'arithmétique des nombres flottants définie par la norme présente quelques particularités mathématiques curieuses et à ne pas oublier :

Il arrive qu'il soit indispensable de réorganiser les calculs pour tenir compte de ces particularités.

Les formats normalisés des flottants

La norme IEEE-754 rev 1985, la plus courante, définit 2 formats principaux de nombres flottants : les nombres en simple précision (le type float du C) et les nombres en double précision (le type double du C). Il existe deux autres types : le simple précision étendue qui est obsolète, et le double précision étendue qui n'est pas encore très utilisé, sauf dans l'ALU des processeurs Intel et autres.

Le type simple précision

Un nombre de type simple précision est codé sur 32 bits en base 2. Sa précision est de 23+1 bits. \(e_{min} \) vaut -126 et \( e_{max} \) vaut +127.

La valeur absolue max en décimal est environ 3,403.1038. Sa précision est d'environ 1,2.10-7.

Le type double précision

Un nombre de type double précision est codé sur 64 bits en base 2. Sa précision est de 52+1 bits. \(e_{min} \) vaut -1022 et \( e_{max} \) vaut +1023.

La valeur absolue max en décimal est environ 1,798.10308. Sa précision est d'environ 2,2.10-16.

Les valeurs spéciales

La représentation de zéro

Il existe deux représentations de zéro dans la norme : +0 et -0. Dans les deux cas, l'exposant et la mantisse sont à zéro. Pour le -0, le bit de signe est à 1. Pour +0, le bit de signe est à 0.

Si vous tentez la comparaison entre +0 et -0, la norme IEEE-754 impose que le résultat soit True.

Les infinis

La norme définit une représentation pour les infinis :

La valeur NaN

Dans l'arithmétique en virgule flottante, certaines opérations sont indéterminées, comme par exemple \( \sqrt{-1} \) ou encore \( \dfrac{0}{0} \) ou encore les opérations sur l'infini.

Dans le cas d'une opération indéterminée, l'arithmétique retourne une valeur spéciale notée NaN, pour Not a Number. Dans la représentation normalisée de NaN, le bit de signe est 0, tous les bits de l'exposant sont à 1, la valeur de la mantisse dépend du compilateur.

Attention à l'algèbre de NaN :

L'arrondi

Définition normalisée

L'arrondi d'un nombre est une opération cruciale en calcul numérique. En effet, le résultat d'une opération entre deux nombres flottants, parfaitement représentables dans l'arithmétique de l'ordinateur est souvent un nombre qui n'est pas représentable dans cette arithmétique. L'exemple le plus simple : 1 et 3 sont représentables sans problème, mais par contre 1/3 ne l'est pas ! Sa mantisse comporte un nombre infini de chiffres. Il faut donc le convertir en un nombre représentable. Il y a deux moyens pour cela : on le tronque à partir d'un certain nombre de décimales, nombre qui dépend du type de la variable qui va stocker le résulat. Ou bien on l'arrondit.

Les propriétés de continuité de la droite des réels nous permet d'affirmer que sur le domaine de définition des nombres flottants n'importe quel réel x est encadré par deux nombres flottants tels que \( x^- \le x \le x^+ \). Si le nombre réel est représentable en flottant, les deux nombres de l'encadrement sont identiques et égaux à la représentation du réel x en flottant.

Par exemple, en simple précision, le rapport 1/3 est encadré par 1.3333333 et 1.3333334. Arrondir un nombre, ici 1/3, cela revient à choisir entre les deux flottants qui l'encadrent. Mais comment choisir et surtout comment faire en sorte que tout le monde choisisse la même règle ?

La norme IEEE-754 définit 4 types d'arrondi et permet de choisir un de ces quatre modes :

La norme IEEE-754 définit la notion d'arrondi "correct" : lorsque le choix de l'arrondi est fait parmi les 4 modes possibles, les 5 opérations aritmétiques courantes, addition, soustraction, multiplication, division et racine carrée donnent un seul résultat possible te toujours reproductible. Ainsi, si le processeur et le compilateur respectent la norme, et si les options de compilation éventuelles permettant de la contourner ne sont pas activées, le programmeur est assuré d'obtenir les mêmes résultats quel que soit l'environnement normalisé. Mais attention, les autres opérations comme les puissances, l'exponentielle, sinus, cos, ne sont pas concernées par la norme. La plupart des compilateurs leur appliquent la norme, mais ce n'est pas obligatoire, alors méfiance...

Même quand le mode d'arrondi est choisi, il faut faire très attention au cumul des arrondis dans les calculs. On parle de propagation des erreurs d'arrondi (comme dans le cas de l'échec de l'interception du SCUD par le Patriot). La norme ne protège pas contre eux et il n'est pas possible de les éviter. Pour les contourner, on en est réduit à borner l'erreur finale, ce qui n'est pas toujours simple.

Quelques exemples instructifs

Sous Python, définissons une constante décimale simple, par exemple 0.15. Et voyons comment Python la retient en mémoire en utilisant la commande format(0.15,'.30f'). J'obtiens le résultat suivant : 0.149999999999999994448884876874. Avant même tout calcul, Python introduit une erreur dans la représentation du nombre décimal 0.15. Ce dernier ne supporte pas une représentation exacte et entraine donc un arrondi.

Autre exemple assez amusant : faites calculer à vos machines l'expression \(12 \left(\dfrac{1}{3} - \dfrac{1}{4} \right) - 1 \), mais faites d'abord le calcul à la main. Il vient très rapidement que le résultat est 0. Python me dit que le résultat est -2.22044604925e-16 . Ma TI89 me dit que le résultat est -4. e-14. Manifestement, les arrondis sur 1/3 et 1/4 n'aident pas nos calculateurs....

Tapez la commande suivante dans votre console Python : format((3.11*2.3)*1.5,'.30f'). Le résultat affiché est 10.729499999999999815258888702374. Maintenant tapez format(3.11*(2.3*1.5),'.30f'), le même calcul mais en déplaçant les parenthèses. Le résultat affiché est : 10.729499999999998038902049302123. Vous constatez avec effarement que la multiplication en arithmétique flottante n'est pas associative ! Même topo pour l'addition. Encore une conséquence de la gestion, pourtant normalisée, des arrondis lors des différentes étapes du calcul dans l'ALU du processeur...

Une dernière chose : dans n'importe quel langage, n'utilisez JAMAIS une condition de sortie de boucle qui serait une comparaison entre un calcul sur des flottants et 0, du style while (x + 2*y) == 0. Vous pouvez être certain du bug annoncé !

La gestion des exceptions

La norme IEEE-754 définit 5 flags qui sont levés lors de la survenue d'une exception : 

Ces flags se propagent dans le calcul et doivent être remis à zéro par le programmeur. Ils sont remis à zéro automatiquement à la fin du process de calcul.

Le programmeur peut accéder à ces flags, par exemple en utilisant les instruction try...except de Python ou les équivalents dans les autres langages. Non seulement il peut, mais il doit ! Chaque calcul qui peut engendrer une exception doit être protégé par un test sur ces flags. Nous avons vu dans le début de cette page ce qui peut résulter d'une division par zéro ou d'un overflow...

Aujourd'hui, tous les langages et tous les compilateurs disposent de moyens pour traiter les exceptions. Les programmeurs n'ont aucune excuse. D'autre part, il est indispensable de vérifier lors des tests de robustesse du code que ces exceptions sont bien traitées. Là aussi, les responsables des tests n'ont aucune excuse : ils disposent de tous les éléments nécessaires.

Des exemples d'erreurs d'arithmétique

Absorption et cancellation

L'arithmétique des nombres flottants normalisée apporte son lot de curiosités, souvent malvenues. Voyons le cas de l'absoprtion et de la cancellation.

Absorption

Le problème se pose quand on tente une addition de deux nombres flottants d'ordre de grandeur très différent. Choisissons par exemple deux flottants simple précision x = 0.1000000e1 et y = 0.5000000e-7 et posons z = x + y.

Essayons de voir ça en pratique, à l'aide de ce petit code en C (compilateur gcc pour Mac OSX sur Intel core i5)

#include <stdio.h>
#include <stdlib.h>

int main(void) {
    float x = 0.1000000e1, y = 0.50000000e-7;

    printf("x = %1.8e\n",x);
    printf("y = %1.8e\n",y);
    printf("x + y = %1.8e\n",x+y);
    return EXIT_SUCCESS;
}

et les résultats obtenus :

x = 1.00000000e+00
y = 5.00000006e-08
x + y = 1.00000000e+00

Au cours du calcul (exactement lorsque le résultat est passé de l'accumulateur de l'ALU vers la mémoire), le nombre y a été absorbé par le nombre x. En choisissant les mêmes valeurs mais en travaillant en double (sur 64 bits) je pallie l'absorption :

x = 1.0000000000000000000e+00
y = 4.9999999999999997737e-08
x + y = 1.0000000499999999182e+00

Notez quand même la valeur que prend la variable y, que l'on croyait avoir défini précisément...

Un autre exemple, prenons deux flottants simple précision :

x = 2.34877563e+02
y = 5.67994416e-01

et calculons la somme x + y, j'obtiens :

x + y = 2.35445557e+02

ce qui n'est pas le résultat recherché !

Cancellation

Une erreur de cancellation (élimination en français, mais l'anglicisme est courant) survient lorsque on soustrait deux nombres flottants de même signe et de valeur très proche.

Utilisons ce petit programme pour calculer l'expression \(\dfrac{1}{x} - \dfrac{1}{x+1}\), avec x = 1.0.106 :

int main(void) {
    float x,y,z;

    x = 1.0e6;
    printf("x = %1.8e\n",x);
    y = (1.0/x) - (1.0/(x + 1));
    printf("y = %1.8e\n",y);

    return EXIT_SUCCESS;
}

J'obtiens le résulat suivant :

y = 9.99999020e-13

ce qui n'est manifestement pas le résultat attendu ! Nous sommes en présence d'une erreur de cancellation. Pour l'éviter, il faut modifier ma formule pour tenir compte des méthodes de calcul de notre processeur et de son arithmétique. Ecrivons la formule équivalente mathématiquement \( \dfrac{x}{x(x+1)}\), sans vouloir simplifier par x. Dans ce cas, j'obtiens :

    z = x/(x*(x+1));
    printf("z = %1.8e\n",z);

avec pour résultat :

z = 9.99998974e-07

ce qui un peu plus conforme à l'attendu ! Une erreur de 6 ordres de grandeur tout de même : imaginez les dégats dans un calcul !

Lorsque vous codez des calculs, ayez présent à l'esprit les ordres de grandeurs des variables et méfiez-vous des valeurs trop différentes pour une addition ou trop proches pour une soustraction.

La suite de Muller

La suite présentée ici est due à Jean-Michel Muller (CNRS/ENS Lyon/INRIA). Elle est souvent utilisée pour illustrer les instabilités numériques propres à l'arithmétique en virgule flottante des ordinateurs.

Son expression est la suivante :

\( \left\{ \begin{array}{lr} u_0 = 2 \\ u_1 = -4 \\ u_n = 111 - \dfrac{1130}{u_{n-1}} + \dfrac{3000}{u_{n-1}u_{n-2}} \\ \end{array} \right. \)

On démontre que la limite théorique de cette suite est 6. Et pourtant, voyons ce que cela donne sur nos ordinateurs. J'ai fait l'essai en C et en Python sur mon Mac. Je vous encourage à essayer sur votre machine, pour voir !

Son traitement en C sur Mac OS

Le code ci-dessus implémente la suite de Muller en C. Sur mon Mac, j'utilise Eclipse comme environnement de développement et le compilateur GNU gcc. Voici le code, on ne peut plus simple :

#include <stdio.h>
#include <stdlib.h>

int main(void) {
    float u0 = 2.0f, u1 = -4.0f;

    printf("u( 0) = %1.19e\n",u0);
    printf("u( 1) = %1.19e\n",u1);

    int n;
    for (n=1;n<= 19;n++)
    {
        float tmp = 111.f - 1130.f/u1 + 3000.f*1/(u1*u0);
        printf("u(%2d) = % 1.19e\n",n+1, tmp);
        u0 = u1;
        u1 = tmp;
    }
    return EXIT_SUCCESS;
}

Notez que je précise explicitement le format float des constantes. Pour rappel, en C le format float est codé sur 4 octets (32 bits) avec une précision de 7 chiffres.

J'obtiens les résultats suivants, pour les vingt premiers termes :

u(0) = 2.0000000000000000000e+00
u(1) = -4.0000000000000000000e+00
u( 2) =  1.8500000000000000000e+01
u( 3) =  9.3783798217773437500e+00
u( 4) =  7.8011646270751953125e+00
u( 5) =  7.1545600891113281250e+00
u( 6) =  6.8088302612304687500e+00
u( 7) =  6.6227531433105468750e+00
u( 8) =  6.9049758911132812500e+00
u( 9) =  1.2952423095703125000e+01
u(10) =  5.7301113128662109375e+01
u(11) =  9.5321716308593750000e+01
u(12) =  9.9694656372070312500e+01
u(13) =  9.9981079101562500000e+01
u(14) =  9.9998840332031250000e+01
u(15) =  9.9999923706054687500e+01
u(16) =  9.9999992370605468750e+01
u(17) =  1.0000000000000000000e+02
u(18) =  1.0000000000000000000e+02
u(19) =  1.0000000000000000000e+02
u(20) =  1.0000000000000000000e+02

Les premiers termes, jusqu'à n = 8, oscillent autour de 6. Puis la suite converge rapidement vers 100. J'ai  essayé en employant des variables de type double, codées sur 64 bits avec 19 chiffres de précision : même résultat.

Son traitement en Python 2.7 sur Mac OS

J'ai transposé le code C ci-dessus en Python :

Indice = []
Valeur = []

# les variables u0 et u1 sont de type float Python ()
u0 = 2.0
u1 = -4.0

print "u( 0) = %21.20f" % u0
print "u( 1) = %21.20f" % u1

Indice.append(0);Valeur.append(u0)
Indice.append(1);Valeur.append(u1)

for i in range(1,19):
    tmp = 111.0 - 1130.0/u1 + 3000.0/(u0*u1)
    print "u(%2d) = %21.20f \n" % (i+1,tmp)
    Indice.append(i)
    Valeur.append(tmp)
    u0 = u1
    u1 = tmp

plt.figure()
plt.grid()
plt.plot(Indice, Valeur)
plt.show()

Voici les résultats obtenus :

u( 0) = 2.00000000000000000000
u( 1) = -4.00000000000000000000
u( 2) = 18.50000000000000000000
u( 3) = 9.37837837837837895449
u( 4) = 7.80115273775216877539
u( 5) = 7.15441448097533339023
u( 6) = 6.80678473692481134094
u( 7) = 6.59263276872179204702
u( 8) = 6.44946593405393286957
u( 9) = 6.34845206074662371520
u(10) = 6.27443866272811590079
u(11) = 6.21869676858216280380
u(12) = 6.17585385581539014765
u(13) = 6.14262717048100626016
u(14) = 6.12024870457015879310
u(15) = 6.16608655959809937031
u(16) = 7.23502116553493124229
u(17) = 22.06207846352579338145
u(18) = 78.57557488787223576310
u(19) = 98.34950312216535905918

La suite diverge pour un indice supérieur, mais elle diverge quand même !

Son traitement en Python 3.6 sur Intel Core i7

Avec le même programme adapté pour Python 3.6 (modification de la fonction print), j'obtiens les résultats suivants :

u( 0) = 2.00000000000000000000
u( 1) = -4.00000000000000000000
u( 2) = 18.50000000000000000000
u( 3) = 9.37837837837837895449
u( 4) = 7.80115273775216877539
u( 5) = 7.15441448097533339023
u( 6) = 6.80678473692481134094
u( 7) = 6.59263276872179204702
u( 8) = 6.44946593405393286957
u( 9) = 6.34845206074662371520
u(10) = 6.27443866272811590079
u(11) = 6.21869676858216280380
u(12) = 6.17585385581539014765
u(13) = 6.14262717048100626016
u(14) = 6.12024870457015879310
u(15) = 6.16608655959809937031
u(16) = 7.23502116553493124229
u(17) = 22.06207846352579338145
u(18) = 78.57557488787223576310
u(19) = 98.34950312216535905918

Les résultats sont identiques à ceux obtenus sur le Mac. Ce n'est guère étonnant, les processeurs et les compilateurs des deux machines respectent la norme IEEE-754.

L'explication

Quelque soit le compilateur et le processeur que vous utiliserez, la limite que vous observerez sera 100. Ce qui montre les limites du calcul sur nos ordinateurs, dont il faut avoir conscience.

Mais pourquoi cette anomalie ? Pour comprendre, il faut réorganiser l'expression du terme de récurrence un comme suit :

\( u_n = \dfrac{\alpha 100^{n+1} + \beta 6^{n+1} + \gamma 5^{n+1}}{\alpha 100^n + \beta 6^n + \gamma 5^n}\)

Les valeurs des paramètres \( \alpha \), \( \beta \) et \( \gamma \) dépendent de la valeur de u0 et u1. Si \( \alpha \) est non nul, la suite converge vers 100, et sinon, avec \( \beta \) non nul, elle converge vers 6.

Dans notre cas avec les valeurs de e u0 et u1, on a  \( \alpha = 0 \), \( \beta = -3  \) et \( \gamma = 4\). \( \alpha \) étant nul, la suite converge mathématiquement vers 6.

Nous avons là un bel exemple d'erreur due à un cumul d'arrondis. Le terme \( \alpha \) devrait être nul tout au long du calcul, mais les erreurs d'arrondis font que \( \alpha \) n'est pas vraiment nul, mais plutôt très petit. A partir d'un certain rang d'indice n, variable selon la précision du calcul, la suite converge vers 100 plutôt que 6.

Ce genre de problème peut vous arriver au détour d'un algorithme. Il est toujours nécessaire de faire une analyse mathématique de nos algorithmes pour en déterminer les éventuelles instabilités numériques et leur sensibilité à l'arithmétique du calculateur.

Un exemple de propagation des erreurs d'arrondi dans un calcul

Soit à calculer numériquement la valeur de e-10, en utilisant la série \(\displaystyle e^{-10} \approx \sum_{k=0}^n (-1)^k \dfrac{10^k}{k!} \).

Pour implémenter cette série et effectuer le calcul, j'utilise le script Python suivant :

from scipy import exp

# définition de la fonction de calcul de la factorielle def fact(n):
    if n<2:
        return 1
    else:
        return n*fact(n-1)
       
# valeur approchée de e-10
print format(exp(-10), '.30f')

# affichage des différents termes de la série
x = 0.0
for k in range(10):
    y = (-1)**k*((10.0**k)/fact(k))
    print 'u(%d) = %.30f' % (k,y,)
    x += y
print u'résultat calculé :%.30f' % (x,)

L'exécution de ce script produit les résultats suivants :

0.000045399929762484854173145571
u(0) = 1.000000000000000000000000000000
u(1) = -10.000000000000000000000000000000
u(2) = 50.000000000000000000000000000000
u(3) = -166.666666666666657192763523198664
u(4) = 416.666666666666685614472953602672
u(5) = -833.333333333333371228945907205343
u(6) = 1388.888888888888914152630604803562
u(7) = -1984.126984126984098111279308795929
u(8) = 2480.158730158730122639099135994911
u(9) = -2755.731922398589176737004891037941
résultat calculé :-1413.144620811287495598662644624710

La première ligne indique la valeur arrondie de e-10. Les lignes suivantes donne la valeur de chaque terme de la série entre 0 et 9. La dernière ligne indique le résultat calculé par la sommation des termes de la série. Le moins que l'on puisse dire, c'est qu'il diffère très sensiblement du résultat attendu ! Il s'agit d'une erreur d'arrondi qui se propage dans le calcul pour donner un résultat aberrant.

Une solution est modifier la méthode de calcul en remarquant que e-10 = 1/e10. Le calcul de e10 utilise aussi une série mais qui ne produit pas d'erreur d'arrondi :

# autre mode de calcul
for k in range(20):
    y = ((10.0**k)/fact(k))
    print 'u(%d) = %.30f' % (k,y,)
    x += y
print u'résultat calculé :%.30f' % (1.0/x,)

Le résultat final devient :
résultat calculé :0.000048692048250898442086128520
ce qui est déjà plus proche du résultat attendu ! Notez que j'ai augmenté le nombre de termes de la série pour améliorer la précision. Si vous faites cela dans le premier calcul, le résultat restera toujours aberrant.

Conclusion

Nous avons la chance de disposer d'une arithmétique en virgule flottante normalisée (IEEE-754) et performante. Elle est certes améliorable (elle l'est continuellement) mais elle permet de programmer avec efficacité.

Toutefois, la conception d'un algorithme de calcul doit être très prudente : il convient de pourchasser toutes les causes possibles d'erreurs, en particulier les erreurs d'arrondi qui se propagent ou pas, les erreurs de troncature, les instabilités numériques, les cancellations et absorptions, et autres...

Il est également nécessaire de traquer tous les dépassements de capacité que les opérations sur des données peuvent engendrer, et respecter la règle essentielle : un code doit traiter les erreurs et ne pas se "planter". L'IEEE-754 le permet : mettez en oeuvre tous les outils à votre disposition.

Toutes ces précautions nécessitent une bonne connaissance de la norme : étudiez-la et faites des manip.

Enfin, comme nul n'est infaillible, il est IMPERATIF de pratiquer tous les tests de robustesse et d'endurance possibles sur vos codes. Ne mettez JAMAIS en service un code, surtout s'il interagit avec des systèmes externes, sans avoir réaliser tous les tests et avoir obtenus des résultats satisfaisants. Pour avoir négligé cet impératif, des catastrophes sont survenues et surviendront encore. Elles sont évitables.

Contenu et design par Dominique Lefebvre - www.tangenteX.com mai 2019 Licence Creative Commons Contact :