Licence CC BY

Aperçu de Gappa, un vérificateur de calculs numériques

Vous savez peut-être que les calculs sur ordinateur peuvent être un peu piégeux. En effet, les approximations des nombres réels utilisées par les ordinateurs réservent des surprises aux non-initiés. Un exemple simple est la comparaison suivante entre nombres à virgule flottante suivant, telle qu’évaluée par l’interpréteur Python :

>>> 0.1 + 0.2 == 0.3
False

Cependant, il existe bien d’autres cas plus subtils.

Ce genre d’entourloupes plus ou moins faciles à identifier, et l’importance cruciale de nombreux codes de calcul a mené des chercheurs à développer un outil pour aider à prouver que les calculs que l’on souhaite faire sont effectivement faits avec la précision attendue. Cet outil s’appelle Gappa et nous donnerons dans cet article un bref aperçu de son utilisation et de son fonctionnement.

Principe de fonctionnement

Utilisation

Il est possible de se procurer Gappa sur le site officiel de l’outil.

Gappa s’utilise en ligne de commande. Il lit en entrée un fichier (ou l’entrée standard) qui contient l’énoncé du problème à résoudre et renvoie le résultat sur la sortie standard, accompagné d’éventuelles erreurs ou avertissements. C’est très simple d’utilisation.

Un script Gappa est structuré généralement comme suit :

  • une définition du problème à résoudre, sous la forme de définitions de paramètres, de calculs, de modes d’arrondis, etc. ;
  • un énoncé de l’objectif à atteindre, souvent sous la forme d’un intervalle à prouver ou à déterminer concernant l’erreur absolue ou des erreurs relatives étant donné des intervalles de départ pour les paramètres du calcul ;
  • des indications pour aider le solveur de Gappa, sous la forme de règles de réécritures, de relations entre des expressions et leurs approximations, ou des indices sur la manière de découper les intervalles, entre autres choses.

Fonctionnement

Gappa fait de l’arithmétique des intervalles. Autrement dit, les valeurs sont représentées par des intervalles qui évoluent au gré des calculs. Gappa détermine rigoureusement l’évolution de ces intervalles dans des suites de calculs, fait de la dichotomie pour éviter des intervalles trop généraux et donc inutiles, mais aussi garde en tête des relations entre expressions pour atteindre une certaine efficacité sans que l’utilisateur ait à guider le travail outre mesure.

Je ne m’étendrais pas plus sur le fonctionnement interne de Gappa, puisque tout est expliqué mieux que je ne saurais le faire dans la bibliographie.

Exemples simples

Un exemple des plus basiques

Demandons à Gappa de vérifier que la somme de 1 et 1 est bien égale à 2. Pour ce faire, je donne à Gappa l’entrée { 1 + 1 = 2 }. Ce qui est entre crochets correspond à la formule logique que Gappa doit tenter de prouver.

Pour l’entrée :

{ 1 + 1 = 2 }

Gappa répond… rien. Cela signifie que tout va bien et que la propriété a été prouvée.

Quand on demande quelque chose que Gappa n’arrive pas à prouver, il nous le fait savoir bruyamment. Par exemple, pour l’entrée :

{ 1 + 1 = 3 }

Gappa répond :

Error: some properties were not satisfied:
  EQL(1 + 1, 3)

Ici, on aperçoit un peu des entrailles de Gappa, qui représente l’égalité sous la forme du prédicat EQL.

Il n’est évidemment pas possible de prouver la précédente égalité, mais Gappa peut également renvoyer une erreur sur des formules vraies s’il ne sait pas les prouver.

Un exemple plus amusant

Là où Gappa devient plus serviable, c’est qu’il peut compléter des intervalles dans des formules logiques de sorte à ce qu’elles soient vraies. Les intervalles ne sont pas forcément optimaux, mais sont corrects à tous les coups.

Par exemple, si on veut savoir dans quel intervalle se trouve x2 pour x compris entre -1 et 2, Gappa peut en trouver un pour nous. L’entrée ici se lit comme « si x est dans l’intervalle [-1, 1], alors x2 est dans quel intervalle ? ».

L’entrée :

{ x in [-1, 2] -> x*x in ? /\ x + 1 in ?}

La réponse :

Results:
  x * x in [0, 4]

On a bien la réponse attendue, à savoir que x2 est compris entre 0 et 4.

Notez bien qu’ici, Gappa considère x comme un réel, vu qu’il n’y a pas de mention spéciale à son sujet.

Prouver la précision d'un algorithme

Position du problème

Suite convergeant vers une racine carrée

L’algorithme auquel nous nous intéressons est une suite qui permet d’approximer la racine carrée d’un nombre. La suite est la suivante :

x0=1x_0 = 1 xn+1=12(xn+axn)x_{n+1} = \frac{1}{2}\left(x_n + \frac{a}{x_n}\right)

aa est le réel dont on cherche la racine carrée.

Il est possible de montrer que la suite converge et que la limite est bien la racine carrée de aa :

limnxn=a\lim_{n \rightarrow \infty} x_n = \sqrt a

Une telle suite ne présente pas de véritable intérêt pratique, parce que bien qu’elle pourrait être utilisée, il existe de meilleurs algorithmes pour extraire des racines carrées.

Nombre d’itérations fini

Pour les besoins de l’exemple, nous nous contenterons de regarder un nombre fini d’itérations, à savoir 3.

Nous voulons savoir quelle précision serait atteignable avec si peu d’itérations pour des nombres compris entre 0,1 et 10.

Un petit calcul numérique préliminaire avec Python me suggère qu’on devrait espérer une erreur absolue comprise entre 0 et environ 0,034 sur cet intervalle.

Première version naïve pour Gappa

En première approche, il est possible de donner à Gappa une formulation du problème la plus simple possible et lui demander de trouver l’intervalle pour l’erreur.

x0 = 1;
x1 = 0.5 * (x0 + a / x0);
x2 = 0.5 * (x1 + a / x1);
x3 = 0.5 * (x2 + a / x2);

{ a in [0.1, 10] -> x3 - sqrt(a) in ?}

La réponse est un peu décevante :

Results:
  x3 - sqrt(a) in [-869305569053949035b-58 {-3.01601, -2^(1.59264)}, 836020376618369155b-55 {23.2042, 2^(4.53632)}]

Gappa nous donne une fenêtre grossière, entre environ -3 et environ 23.

Deuxième version, avec des indications de résolution

Il est possible de donner des indications à Gappa pour l’aider dans sa recherche d’intervalle.

x0 = 1;
x1 = 0.5 * (x0 + a / x0);
x2 = 0.5 * (x1 + a / x1);
x3 = 0.5 * (x2 + a / x2);

{ a in [0.1, 10] -> x3 - sqrt(a) in ?}

$ a in 100;

La dernière ligne $ a in 100; indique à Gappa qu’il devrait effectuer sa recherche en découpant l’intervalle dans lequel se trouve aa en tranches, ici en 100 parties. C’est une heuristique qui peut aider, parce que l’algorithme inverse des valeurs inférieures à 1, qui vont avoir tendance à créer des nombres relativement grands et donc élargir excessivement les bornes estimées par Gappa.

La réponse est plus satisfaisante :

Results:
  x3 - sqrt(a) in [-375550141465633147b-61 {-0.162869, -2^(-2.61822)}, 465140072240540851b-61 {0.201722, 2^(-2.30956)}]

Ici, l’erreur est estimée entre environ -0,16 et environ 0,2. C’est assez loin de ce qu’on peut attendre encore.

On peut augmenter le découpage. Si on donne 10 000, l’intervalle est réduit à environ -0,001 et 0,034. La borne inférieure n’est pas incroyable, vu qu’on peut attendre à quelque chose de très proche de zéro. La borne supérieure est proche du meilleur qu’on peut avoir, mais avoir à découper autant n’est pas très satisfaisant.

Il est possible de donner d’autres types d’indices à Gappa

Troisième version, avec des règles de réécriture

Il est possible d’aider Gappa en l’aidant à identifier les termes d’erreur, à savoir les choses de la forme xax - \sqrt a.

Cela se fait en utilisant des règles de réécriture. Ici, on indique la chose suivante pour chaque itération :

12(xn+axn)a=(xna)22xn\frac{1}{2}\left(x_n + \frac{a}{x_n}\right) - \sqrt a = \frac{\left(x_n -\sqrt a \right)^2}{2\cdot x_n}

Cela peut aider Gappa à apparenter les termes et mieux gérer son inférence. Cette idée est directement inspirée d’un exemple de la documentation de Gappa.

x0 = 1;
x1 = 0.5 * (x0 + a / x0);
x2 = 0.5 * (x1 + a / x1);
x3 = 0.5 * (x2 + a / x2);

{ a in [0.1, 10] -> x3 - sqrt(a) in ?}

0.5 * (x0 + a / x0) - sqrt(a) -> 0.5 / x0 * (x0 - sqrt(a)) * (x0 - sqrt(a));
0.5 * (x1 + a / x1) - sqrt(a) -> 0.5 / x1 * (x1 - sqrt(a)) * (x1 - sqrt(a));
0.5 * (x2 + a / x2) - sqrt(a) -> 0.5 / x2 * (x2 - sqrt(a)) * (x2 - sqrt(a));
$ a in 100;

Demander de découper l’intervalle de a en tranches fines est encore nécessaire pour s’approcher de la bonne réponse, mais avec seulement 100, on obtient un résultat proche de ce qu’on avait avec 10 000 tout à l’heure, alors qu’on découpe quand même 100 fois moins. Ici, on se trouve avec une erreur comprise entre -10-14 (presque zéro, donc) et 0,034 et des poussières, ce qui est satisfaisant par rapport à nos attentes.

Warning: x1 - sqrt(a) and 5e-1 / x0 * (x0 - sqrt(a)) * (x0 - sqrt(a)) are not trivially equal.
         Difference: ((a) - (sqrt(a))^2) / (2)
Warning: the expression (a) + 1 has been assumed to be nonzero when checking a rewriting rule.
Warning: the expression (a) + 1 has been assumed to be nonzero when checking a rewriting rule.
Warning: x2 - sqrt(a) and 5e-1 / x1 * (x1 - sqrt(a)) * (x1 - sqrt(a)) are not trivially equal.
         Difference: ((a) + (a)^2 - (a) * (sqrt(a))^2 - (sqrt(a))^2) / (2 * (a) + 1 + (a)^2)
Warning: the expression (a) + 1 has been assumed to be nonzero when checking a rewriting rule.
Warning: the expression 6 * (a) + (a)^2 + 1 has been assumed to be nonzero when checking a rewriting rule.
Warning: the expression 6 * (a) + (a)^2 + 1 has been assumed to be nonzero when checking a rewriting rule.
Warning: x3 - sqrt(a) and 5e-1 / x2 * (x2 - sqrt(a)) * (x2 - sqrt(a)) are not trivially equal.
         Difference: (-2 * (a)^6 * (sqrt(a))^2 - 62 * (a)^2 * (sqrt(a))^2 - 88 * (a)^3 * (sqrt(a))^2 + 2 * (a) + 20 * (a)^2 + 62 * (a)^3 - 2 * (sqrt(a))^2 + 88 * (a)^4 - 62 * (a)^4 * (sqrt(a))^2 + 62 * (a)^5 + 20 * (a)^6 + 2 * (a)^7 - 20 * (a)^5 * (sqrt(a))^2 - 20 * (a) * (sqrt(a))^2) / (15 * (a) + 1 + 77 * (a)^2 + 163 * (a)^3 + 163 * (a)^4 + 77 * (a)^5 + 15 * (a)^6 + (a)^7)
Results:
  x3 - sqrt(a) in [-930798660757027715b-106 {-1.1473e-14, -2^(-46.3088)}, 638742597432983727b-64 {0.0346263, 2^(-4.85199)}]

Gappa indique plein d’avertissements ! Certaines vérifications sont faites par Gappa pour éviter des erreurs grossières de réécriture, mais il ne sait pas prouver la correction de toutes les réécritures, surtout s’il y a besoin d’hypothèses non-indiquées pour le faire. Dans notre cas, il ne peut pas s’assurer que a(a)2a - (\sqrt a)^2 est égal à zéro par exemple, ou que certaines valeurs apparaissant dans des divisions sont effectivement non-nulles.

Et avec des flottants ?

Gappa connaît bien les flottants. Ici, je lui ai indiqué de considérer ar comme un réel et a comme un flottant, mais surtout d’effectuer toutes les opérations en utilisant des flottants binary64 avec un arrondi au pair le plus proche, autrement dit, les flottants usuels de mon ordinateur.

@round = float<ieee_64,ne>;
a = round(ar);
x0 round= 1;
x1 round= 0.5 * (x0 + a / x0);
x2 round= 0.5 * (x1 + a / x1);
x3 round= 0.5 * (x2 + a / x2);

{ ar in [0.1, 10] -> x3 - sqrt(ar) in ?}

0.5 * (x0 + a / x0) - sqrt(a) -> 0.5 / x0 * (x0 - sqrt(a)) * (x0 - sqrt(a));
0.5 * (x1 + a / x1) - sqrt(a) -> 0.5 / x1 * (x1 - sqrt(a)) * (x1 - sqrt(a));
0.5 * (x2 + a / x2) - sqrt(a) -> 0.5 / x2 * (x2 - sqrt(a)) * (x2 - sqrt(a));
$ a in 100;

Tout se passe bien pour nous, parce qu’après les avertissements, on a toujours un résultat satisfaisant par rapport à l’étape d’avant :

[...]
Results:
  x3 - sqrt(aa) in [-236926282298699135b-104 {-1.16814e-14, -2^(-46.2828)}, 638742597432994563b-64 {0.0346263, 2^(-4.85199)}]

Voilà ! Ce n’était qu’un petit aperçu de cet outil. Bien que de prise en main simple, il peut s’avérer un peu complexe à dompter, afin d’en tirer le maximum.

Au-delà d’une utilisation comme outil indépendant, Gappa peut s’interfacer avec d’autres outils de preuve formelle. Par exemple, Gappa peut générer des preuves destinées à être vérifiées par le moteur de Coq, ou encore s’interfacer avec la plateforme Why (et donc Frama-C par extension, entre autres).

Références

5 commentaires

Tiens, je reconnais des noms dans la liste bibliographique de cet outil. Il faut vraiment que je me mette à Coq moi !

+0 -0

En outil en cours de développement pour les preuves sur les flottants, un de mes collègues bosse sur Colibri2 qui est dispo via Opam depuis quelques mois.

Tiens, je reconnais des noms dans la liste bibliographique de cet outil. Il faut vraiment que je me mette à Coq moi !

ache

My 2 cents : si tu veux prouver des algorithmes, attaque plutôt par Why3.

My 2 cents : si tu veux prouver des algorithmes, attaque plutôt par Why3.

En vrai, je me demande pourquoi (sans jeu de mot). ^^

+1 -0

Je pense que ce n’est pas une bonne idée que je fasse un tel billet. Je ne serais pas du tout impartial :D .

Mais dans l’idée, le principal se résume au fait que c’est un langage de prog et de spécification qui a été conçu dès le départ pour prouver des programmes, et pas des théorèmes. La spécification n’est par ailleurs pas transmise comme des types dépendants, ce qui donne à l’usage un peu plus de flexibilité sur ce que tu spécifies ou pas et ce que tu prouves ou pas. Conséquence au niveau de la preuve, on ne cherche pas à expliquer toutes les étapes du typages, donc on a accès à une pelleté de SMT solvers disponibles qui évitent de faire des tas d’étapes de preuves pas forcément très intéressantes. Comme il y a du code ghost, tu peux quand même aider la preuve de manière semi interactive. Et si vraiment tu n’arrives pas à faire ta preuve comme ça, tu peux quand même faire du Coq pour terminer ta preuve.

Connectez-vous pour pouvoir poster un message.
Connexion

Pas encore membre ?

Créez un compte en une minute pour profiter pleinement de toutes les fonctionnalités de Zeste de Savoir. Ici, tout est gratuit et sans publicité.
Créer un compte