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

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.

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 x² 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, 2], alors x² est dans quel intervalle ? ».

L’entrée :

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

La réponse :

Results:
  x * x in [0, 4]

On a bien la réponse attendue, à savoir que x² 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.

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 :

$x_0 = 1$ $x_{n+1} = \frac{1}{2}\left(x_n + \frac{a}{x_n}\right)$

où $a$ 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 $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 $a$ 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 s’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 $x - \sqrt a$.

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

$\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 - (\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)}]