J'ai été surpris de constater qu'il est difficile de trouver des informations sur ce problème en russe, comme si peu de gens se soucient du fait que les bibliothèques mathématiques utilisées dans les compilateurs modernes ne donnent parfois pas un résultat correctement arrondi. Je m'inquiète de cette situation, car je ne travaille que sur le développement de telles bibliothèques mathématiques. Dans la littérature étrangère, ce problème est bien couvert, j'ai donc décidé de le présenter en russe sous une forme scientifique populaire, en m'appuyant sur des sources occidentales et encore un peu d'expérience personnelle.

---

Amis, pour votre commodité, l'article est également disponible sous forme de présentation vidéo (près de 34 minutes), ce format est plus adapté aux lecteurs qui ont du mal à construire les images mathématiques nécessaires dans leur tête, car il y a beaucoup de matériel illustratif dans la présentation. Les informations sur la vidéo sont complètement identiques au contenu de l'article. Veuillez agir à votre convenance.

---

Je répète qu'il ne s'agit pas d'un article scientifique, mais d'un article de vulgarisation scientifique, après l'avoir lu, vous vous familiariserez brièvement avec cela.

• Les fonctions élémentaires transcendantales (exp, sin, log, cosh et autres) qui fonctionnent avec l'arithmétique en virgule flottante sont arrondies incorrectement, parfois elles font une erreur dans le dernier bit.
• La raison des erreurs ne réside pas toujours dans la paresse ou les faibles qualifications des développeurs, mais dans une circonstance fondamentale, que la science moderne n'a pas encore été en mesure de surmonter.
• «», - .
• , , , , exp2(x) pow(2.0, x).

Pour comprendre le contenu de cet article, vous devez vous familiariser avec le format à virgule flottante IEEE-754. Ce sera suffisant si vous comprenez au moins simplement que, par exemple, il s'agit de: 0x400921FB54442D18 - numéro pi au format double précision (binaire64 ou double), c'est-à-dire que vous comprenez simplement ce que je veux dire par cet enregistrement; Je n'exige pas de pouvoir effectuer de telles transformations à la volée. Et je vais vous rappeler les modes d'arrondi dans cet article, c'est une partie importante de l'histoire. Il est également souhaitable de connaître l'anglais "programmeur", car il y aura des termes et des citations de la littérature occidentale, mais vous pouvez vous en tirer avec un traducteur en ligne.



Exemples d'abord, afin que vous compreniez immédiatement quel est le sujet de la conversation. Maintenant, je vais donner le code en C ++, mais si ce n'est pas votre langage, alors je suis sûr que vous comprendrez toujours facilement ce qui est écrit. Veuillez regarder ce code:

#include <stdio.h>
#include <cmath>

int main() {
  float x = 0.00296957581304013729095458984375f;  // ,  .
  float z;
  z = exp2f(x);  // z = 2**x  .
  printf ("%.8f\n", z);  //      8   .
  z = powf(2.0f, x);  // z = 2**x  
  printf ("%.8f\n", z);  //   .
  return 0;
}

Le nombre x est écrit intentionnellement avec un tel nombre de chiffres significatifs afin qu'il soit exactement représentable dans le type float, c'est-à-dire que le compilateur le convertit en code binaire sans arrondi. Après tout, vous savez très bien que certains compilateurs ne sont pas capables d'arrondir sans erreurs (si vous ne savez pas, indiquez dans les commentaires, j'écrirai un article séparé avec des exemples). Ensuite dans le programme, nous devons calculer 2 ^x , mais faisons-le de deux façons: la fonction exp2f (x) et l'exponentiation explicite de deux powf (2.0f, x). Le résultat, bien sûr, sera différent, car j'ai dit plus haut que les fonctions élémentaires ne peuvent pas fonctionner correctement dans tous les cas, et j'ai spécialement sélectionné un exemple pour le montrer. Voici le résultat:

1.00206053
1.00206041

Quatre compilateurs m'ont donné ces valeurs: Microsoft C ++ (19.00.23026), Intel C ++ 15.0, GCC (6.3.0) et Clang (3.7.0). Ils diffèrent par un bit le moins significatif. Voici le code hexadécimal de ces nombres:

0x3F804385  // 
0x3F804384  // 

Souvenez-vous de cet exemple, c'est là-dessus que nous nous pencherons sur l'essence du problème un peu plus tard, mais pour l'instant, afin que vous ayez une impression plus claire, veuillez consulter les exemples pour le type de données double précision (double, binary64) avec d'autres fonctions élémentaires. Je présente les résultats dans le tableau. Les réponses correctes (lorsqu'elles sont disponibles) ont * à la fin.







J'espère que vous n'avez pas l'impression que j'ai délibérément fait des tests complètement uniques que vous pouvez à peine trouver? Si tel est le cas, cuisinons à genoux une énumération complète de tous les arguments fractionnaires possibles pour la fonction 2 ^x pour le type de données float. Il est clair que nous ne nous intéressons qu'aux valeurs de x comprises entre 0 et 1, car d'autres arguments donneront un résultat qui ne diffère que par la valeur dans le champ d'exposant et ne sont pas intéressants. Vous comprenez vous-même:

$$

$$2^x = 2^{[x]}\cdot2^{\{x\}}.$$

Après avoir écrit un tel programme (le texte caché sera ci-dessous), j'ai vérifié la fonction exp2f et le nombre de valeurs erronées qu'elle produit sur l'intervalle x de 0 à 1.







Il ressort clairement du programme ci-dessous que le nombre d'arguments x testés était 197612997. Il s'avère que, par exemple, Microsoft C ++ calcule incorrectement la fonction 2 ^x pour près d'un pour cent d'entre eux. Ne vous réjouissez pas, chers fans de GCC et Clang, c'est juste que cette fonction est correctement implémentée dans ces compilateurs, mais pleine d'erreurs dans d'autres.





Je n'ennuierai pas le lecteur avec ces exemples, l'essentiel ici était de montrer que les implémentations modernes de fonctions transcendantales peuvent arrondir le dernier bit de manière incorrecte, et que différents compilateurs font des erreurs à différents endroits, mais aucun d'entre eux ne fonctionnera correctement. Au fait, la norme IEEE-754 autorise cette erreur dans le dernier bit (dont je parlerai ci-dessous), mais cela me semble toujours étrange: ok double, c'est un type de données volumineux, mais float peut être vérifié par force brute! Était-ce si difficile à faire? Pas difficile du tout, et j'ai déjà montré un exemple.



Notre code d'énumération contient une fonction "auto-écrite" de calcul correct 2 ^xen utilisant un polynôme approximatif du 10e degré, et il a été écrit en quelques minutes, car ces polynômes sont dérivés automatiquement, par exemple, dans le système d'algèbre informatique de Maple. Il suffit de définir une condition pour que le polynôme fournisse 54 bits de précision (pour cette fonction, 2 ^x ). Pourquoi 54? Mais vous le découvrirez bientôt, juste après que je vous ai dit l'essence du problème et que je vous dise pourquoi, en principe, il est maintenant impossible de créer des fonctions transcendantales rapides et correctes pour le type de données de quadruple précision (binaire128), bien qu'il y ait déjà des tentatives pour attaquer ce problème en théorie.

Arrondi par défaut et son problème

Si vous n'êtes pas plongé dans le développement de bibliothèques mathématiques, alors il n'y a rien de mal à oublier la règle d'arrondi par défaut pour les nombres à virgule flottante selon la norme IEEE-754. Par conséquent, je vais vous le rappeler. Si vous vous souvenez bien de tout, regardez au moins la fin de cette section de toute façon, vous êtes surpris: je vais vous montrer une situation où arrondir un nombre peut être très difficile.



Vous pouvez facilement vous rappeler ce que signifie «arrondir» (à plus l'infini), «arrondir» (à moins l'infini) ou «arrondir à zéro» par le nom (le cas échéant, il y a Wikipedia). Les principales difficultés pour les programmeurs sont d'arrondir "au plus proche, mais dans le cas d'une distance égale du plus proche - à celui avec le dernier chiffre pair". Oui, c'est ainsi que se traduit ce mode d'arrondi, que la littérature occidentale appelle en bref: «Arrondir les liens les plus proches à pair».



Ce mode d'arrondi est utilisé par défaut et fonctionne comme suit. Si, à la suite de calculs, la longueur de la mantisse s'est avérée supérieure à ce que le type de données résultant peut accepter, l'arrondi est effectué à la plus proche de deux valeurs possibles. Cependant, une situation peut survenir lorsque le nombre d'origine s'est avéré être exactement au milieu entre les deux plus proches, alors le résultat est choisi pour lequel le dernier bit (après arrondi) s'avère pair, c'est-à-dire égal à zéro. Prenons quatre exemples où vous devez arrondir à deux bits après le point décimal binaire:

1. Arrondir 1,00 1 001. Le troisième bit après la virgule décimale est 1, mais il y a un autre 6e bit, qui est 1, ce qui signifie que l'arrondi sera supérieur, car le nombre d'origine est plus proche de 1,01 que de 1,00.
2. 1,001000. , 1,00 1,01, .
3. 1,011000. 1,01 1,10. , .
4. 1,010111. , 1,01, 1,10.

D'après ces exemples, il peut sembler que tout est simple, mais ce n'est pas le cas. Le fait est que parfois nous ne pouvons pas dire avec certitude si nous sommes vraiment au milieu entre deux valeurs. Voir un exemple. Essayons à nouveau d'arrondir à deux bits après la virgule décimale:



1.00 1 0000000000000000000000000000000000001



Il est maintenant évident pour vous que l'arrondi doit être supérieur, c'est-à-dire au nombre 1.01. Cependant, vous regardez un nombre avec 40 bits après la virgule décimale. Et si votre algorithme ne pouvait pas fournir 40 bits de précision et n'atteindrait que 30 bits? Ensuite, il donnera un autre numéro:



1.00 1 000000000000000000000000000



Ne sachant pas qu'à la 40e position (que l'algorithme est incapable de calculer), il y en aura une chérie, vous arrondissez ce nombre à la baisse et obtenez 1,00, ce qui est faux. Vous avez mal arrondi la dernière partie - c'est le sujet de notre discussion. D'après ce qui précède, il s'avère que pour obtenir uniquement le 2ème bit correct, vous devez calculer la fonction jusqu'à 40 bits! Hou la la! Et si la "locomotive" des zéros s'avère encore plus longue? C'est ce dont nous parlerons dans la section suivante.



En passant, c'est l'erreur de nombreux compilateurs lors de la conversion de la notation décimale d'un nombre à virgule flottante au format binaire résultant. Si le nombre décimal d'origine dans le code du programme est trop proche du milieu entre deux valeurs binaires représentables avec précision, il ne sera pas arrondi correctement. Mais ce n'est pas le sujet de cet article, mais une raison pour une histoire séparée.

L'essence du problème d'arrondir le dernier bit significatif

Le problème se manifeste pour deux raisons. Le premier est le rejet délibéré des calculs chronophages en faveur de la vitesse. Dans ce cas, tant que la précision spécifiée est observée, et quels bits il y aura dans la réponse est une question secondaire. La deuxième raison est le dilemme de Table Maker, qui est le sujet principal de notre conversation. Examinons les deux raisons plus en détail.

Première raison

Vous comprenez bien sûr que le calcul des fonctions transcendantales est implémenté par certaines méthodes approximatives, par exemple, par la méthode d'approximation des polynômes ou même (rarement) par développement en série. Pour faire les calculs le plus rapidement possible, les développeurs conviennent d'effectuer le moins d'itérations de la méthode numérique possible (ou de prendre un polynôme du moindre degré possible), tant que l'algorithme autorise une erreur ne dépassant pas la moitié de la valeur du dernier bit de la mantisse. Dans la littérature, cela s'écrit 0.5ulp (ulp = unité en dernier lieu ).



Par exemple, si nous parlons d'un nombre x de type float dans l'intervalle (0,5; 1), la valeur ulp = 2 ^-23 . Sur l'intervalle (1; 2) ulp = 2 ^-22 . En d'autres termes, si x est sur l'intervalle (0; 1) alors 2 ^xsera sur l'intervalle (1,2), et pour assurer une précision de 0.5ulp, il faut, grosso modo, choisir EPS = 2 ^-23 (comme nous désignerons la constante "epsilon", chez les gens ordinaires appelés "erreur", ou "précision", à qui comme vous le souhaitez, ne trouvez pas de faute, s'il vous plaît).



Pour les calculs appliqués, cela suffit, mais le fait que les derniers bits puissent ne pas coïncider avec le résultat absolu n'est pas important pour presque 100% des programmeurs, car il est important pour eux non pas quels seront les bits, mais quelle sera la précision.



Pour ceux qui ne comprennent pas, je vais donner un exemple dans le système de nombres décimaux. Voici deux nombres: 1.999999 et 2.0. Disons que le premier est ce que le programmeur a reçu, et le second est la norme de ce qui aurait dû être obtenu si nous avions des possibilités illimitées. La différence entre eux n'est que d'un millionième, c'est-à-dire que la réponse est calculée avec une erreur de EPS = 10 ^-6 . Cependant, il n'y a pas un seul chiffre correct dans cette réponse. Est-il mauvais? Non, du point de vue du programme d'application, c'est violet, le programmeur arrondira la réponse, disons, à deux décimales et recevra 2,00 (par exemple, il s'agissait de monnaie, 2,00 $), il n'a pas besoin de plus, mais le fait qu'il mettre EPS = 10 ^-6 dans mon programme , puis bien fait, a pris une marge pour l'erreur de calculs intermédiaires et a résolu le problème correctement.



En d'autres termes, ne soyez pas confus: la précision et le nombre de bits (ou chiffres) corrects sont deux choses différentes. Ceux qui ont besoin de précision (c'est presque 100% des programmeurs), le problème discuté ne les concerne pas du tout. Quiconque a besoin d'une séquence de bits pour correspondre à un standard correctement arrondi s'inquiète beaucoup de ce problème, par exemple, les développeurs de bibliothèques de fonctions élémentaires. Néanmoins, il est utile que tout le monde le sache pour le développement général.



Permettez-moi de vous rappeler que c'était la première direction du problème: les derniers éléments de la réponse peuvent être faux car il s'agit d'une solution intentionnelle. L'essentiel est de conserver une précision de 0,5ulp (ou plus). Par conséquent, l'algorithme numérique est sélectionné uniquement à partir de cette condition, si seulement il fonctionne extrêmement rapidement. Dans le même temps, la norme permet la mise en œuvre de fonctions élémentaires sans arrondi correct du dernier bit. Je cite [1, section 12.1] (anglais):

La version 1985 de la norme IEEE 754 pour l'arithmétique à virgule flottante ne spécifiait rien concernant la fonction élémentaire. En effet, on a cru pendant des années que des fonctions correctement arrondies seraient beaucoup trop lentes au moins pour certains arguments d'entrée. La situation a changé depuis et la version 2008 de la norme recommande (mais n'exige pas) que certaines fonctions soient correctement arrondies.

Les fonctions suivantes sont recommandées mais n'ont pas besoin d'être arrondies correctement:

La deuxième raison

Enfin, nous sommes arrivés au sujet de la conversation: le dilemme de Table Maker (abrégé en TMD). Je n'ai pas pu traduire correctement son nom en russe, il a été introduit par William Kahan (père fondateur de l'IEEE-754) dans l'article [2]. Peut-être que si vous lisez l'article, vous comprendrez pourquoi le nom est exactement cela. En bref, l'essence du dilemme est que nous devons obtenir un arrondi absolument précis de la fonction z = f (x), comme si nous avions un enregistrement de bits infini du résultat parfaitement calculé z à notre disposition. Mais il est clair pour tout le monde que nous ne pouvons pas obtenir une séquence infinie. Combien de bits prendre alors? Ci-dessus, j'ai montré un exemple où nous avons besoin de voir 40 bits du résultat afin d'obtenir au moins 2 bits corrects après arrondi. Et l'essence du problème TMD est que nous ne savons pas à l'avance, jusqu'à combien de bits pour calculer la valeur de z afin d'obtenir autant de bits corrects après arrondi que nécessaire. Et s'il y en avait cent ou mille? On ne sait pas à l'avance!



Par exemple, comme je l'ai dit, pour la fonction 2 ^x , pour le type de données float, où la partie fractionnaire de la mantisse n'a que 23 bits, nous devons effectuer le calcul avec une précision de 2 à ⁵⁴ afin que l'arrondi se produise correctement pour tous les arguments x possibles sans exception. Il n'est pas difficile d'obtenir cette estimation par une recherche exhaustive, mais pour la plupart des autres fonctions, en particulier pour les types double ou long double (mettez "classe" si vous savez ce que c'est), ces estimations sont inconnues .



Comprenons déjà pourquoi cela se produit. J'ai délibérément donné le tout premier exemple de cet article avec le type de données float et vous ai demandé de vous en souvenir, car dans ce type, il n'y a que 32 bits et il sera plus facile de le regarder, dans d'autres types de données, la situation est similaire.



Nous avons commencé avec le nombre x = 0,00296957581304013729095458984375, il s'agit d'un nombre exactement représentable dans le type de données float, c'est-à-dire qu'il est écrit de manière à pouvoir être converti en système flottant binaire sans arrondi. Nous calculons 2 ^x , et si nous avions une calculatrice avec une précision infinie, alors nous devrions obtenir (pour que vous puissiez me vérifier, le calcul se fait dans le système WolframAlpha en ligne ):



1.0020604729652405753669743044108123031635398201893943954577320057 ...



Traduisons ce nombre en binaire, disons que 64 bits suffiront:



1.00000000100001110000100 1 000000000000000000000000000001101111101



Le bit d'arrondi (24e bit après la virgule décimale) est souligné. Question: où arrondir? Haut ou bas? Clairement, en haut, vous le savez parce que vous voyez suffisamment de bits et que vous pouvez prendre une décision. Mais regardez attentivement ...



Après le bit d'arrondi, nous avons 29 zéros. Cela signifie que nous sommes très, très proches du milieu entre les deux nombres les plus proches et qu'il suffit de descendre un peu, car la direction de l'arrondi changera. Mais la question est: où sera ce changement? L'algorithme numérique peut séquentiellement, pas à pas, approcher la valeur exacte de différents côtés, et jusqu'à ce que nous passions tous ces 29 zéros et atteignions une précision qui dépasse la valeur du tout dernier zéro dans cette "locomotive", nous ne connaîtrons pas le sens d'arrondi ... Et si, en fait, la bonne réponse devait être:



1.00000000100001110000100 0 11111111111111111111111111111?



Ensuite, l'arrondi sera à la baisse.



Nous ne le savons pas jusqu'à ce que notre précision atteigne le 54e bit après la virgule décimale. Lorsque le 54e bit sera connu exactement, nous saurons exactement de quel nombre des deux nombres les plus proches nous sommes réellement plus proches. Ces nombres sont appelés les points les plus difficiles à arrondir [1, section 12.3] (points critiques pour l'arrondissement), et le nombre 54 est appelé dureté à arrondir, et est désigné par la lettre m dans le livre cité.



La complexité de l'arrondi (m) est le nombre de bits qui est le minimum nécessaire pour garantir que pour tous les arguments d'une certaine fonction f (x) et pour une plage présélectionnée, la fonction f (x) est arrondie correctement au dernier bit (pour différents modes d'arrondi, il peut y avoir différents valeur m). En d'autres termes, pour le type de données float et pour l'argument x de la plage (0; 1) pour le mode d'arrondi "le plus proche pair", la complexité d'arrondi est m = 54. Cela signifie que pour absolument tous les x de l'intervalle (0; 1), nous pouvons mettre dans l'algorithme la même précision ESP = ^2-54 , et tous les résultats seront arrondis correctement à 23 bits après la virgule décimale binaire.



En fait, certains algorithmes sont capables de fournir un résultat exact et basé sur 53 et même 52 bits, la force brute le montre, mais théoriquement, il vous en faut exactement 54. S'il n'y avait pas la possibilité de faire sortir la force brute, nous ne pourrions pas "tricher" et économisez quelques bits, comme je l'ai fait dans le programme de force brute ci-dessus. J'ai pris un polynôme avec un degré inférieur à ce qu'il devrait, mais ça marche toujours, juste parce que j'ai eu de la chance.



Ainsi, quel que soit le mode d'arrondi, nous avons deux situations possibles: soit une «locomotive à vapeur» de zéros surgit dans la zone d'arrondi, soit une «locomotive à vapeur» de uns. La tâche de l'algorithme correct pour le calcul de la fonction transcendantale f (x) est d'affiner la valeur de cette fonction jusqu'à ce que la précision dépasse la valeur du dernier bit de cette "locomotive", et jusqu'à ce qu'il devienne clair qu'à la suite des fluctuations ultérieures de l'algorithme numérique pour le calcul de f (x) les zéros ne se transforment pas en uns, ou vice versa. Dès que tout s'est stabilisé et que l'algorithme a atteint une telle précision qui dépasse les limites de la «locomotive à vapeur», on peut arrondir comme si l'on avait un nombre infini de bits. Et cet arrondi sera avec le dernier bit correct. Mais comment y parvenir?

"Béquilles"

Comme mentionné, le principal problème est d'obtenir l'algorithme pour surmonter la locomotive des zéros ou des uns qui vient immédiatement après le bit d'arrondi. Lorsque la locomotive est vaincue et que nous la voyons dans son ensemble, cela équivaut au fait que ces zéros ou ces uns ont déjà été calculés exactement , et nous savons déjà exactement dans quelle direction l'arrondi se produira maintenant. Mais si nous ne connaissons pas la longueur de la locomotive, comment pouvons-nous concevoir un algorithme?

La première "béquille"

Il peut sembler au lecteur que la réponse est évidente: prenez l'arithmétique avec une précision infinie et mettez un nombre délibérément excessif de bits, et si cela ne suffit pas, mettez-en une autre et recalculez. En général, c'est correct. Ceci est fait lorsque la vitesse et les ressources de l'ordinateur ne jouent pas un rôle particulier. Cette approche a un nom: la stratégie multiniveau de Ziv [1, section 12.3]. Son essence est extrêmement simple. L'algorithme doit prendre en charge les calculs à plusieurs niveaux: un calcul préliminaire rapide (dans la plupart des cas, il s'avère être final), un calcul plus lent mais plus précis (enregistre dans la plupart des cas critiques), encore plus lent, mais encore plus précis (quand c'est absolument "mauvais "Je devais le faire) et ainsi de suite.



Dans la très grande majorité des cas, il suffit de prendre la précision légèrement supérieure à 0,5ulp, mais si une "locomotive" apparaît, alors on l'augmente. Tant que la «locomotive à vapeur» demeure, nous augmentons la précision jusqu'à ce qu'il soit strictement clair que de nouvelles fluctuations de la méthode numérique n'affecteront pas cette «locomotive à vapeur». Ainsi, par exemple, dans notre cas, si nous avons atteint ESP = ^2-54 , alors à la 54ème position apparaît une unité qui, en quelque sorte, "protège" la locomotive des zéros et garantit qu'il n'y aura plus de soustraction d'une valeur supérieure ou égale à ^2-53 et les zéros ne deviendront pas des uns, faisant glisser le bit d'arrondi à zéro avec lui.



C'était une présentation scientifique populaire, tout de même avec le test d'arrondi de Ziv, où il est montré à quelle vitesse, en une seule étape, pour vérifier si nous avons atteint la précision souhaitée, peut être lu dans [1, Chapitre 12], ou dans [3, Section 10,5].



Le problème avec cette approche est évident. Il est nécessaire de concevoir un algorithme de calcul de chaque fonction transcendantale f (x) afin qu'au cours de la pièce, il soit possible d'augmenter la précision des calculs. Pour la mise en œuvre logicielle, ce n'est toujours pas si effrayant, par exemple, la méthode de Newton permet, grosso modo, de doubler le nombre de bits exacts après la virgule décimale à chaque itération. Vous pouvez doubler jusqu'à ce que cela devienne "assez", bien que ce soit un processus assez long, je dois admettre que la méthode de Newton n'est pas toujours justifiée, car elle nécessite de calculer la fonction inverse f ^-1(x), qui dans certains cas peut ne pas être plus simple que de calculer f (x) lui-même. Pour la mise en œuvre matérielle, la "stratégie Ziva" est totalement inadaptée. L'algorithme, câblé dans le processeur, doit effectuer une série d'actions avec le nombre de bits déjà prédéfini, et c'est assez problématique à mettre en œuvre si on ne connaît pas ce nombre à l'avance. Faire l'inventaire? Combien?



L'approche probabiliste pour résoudre le problème [1, section 12.6] nous permet d'estimer la valeur de m (rappelez-vous, c'est le nombre de bits, ce qui est suffisant pour un arrondi correct). Il s'avère que la longueur de la "locomotive" au sens probabiliste est légèrement supérieure à la longueur de la mantisse du nombre. Ainsi, dans la plupart des cas, il suffira de prendre m un peu plus du double de la valeur de la mantisse, et ce n'est que dans de très rares cas qu'il sera nécessaire d'en prendre encore plus. Je cite les auteurs de ce travail: "on en déduit qu'en pratique, m doit être légèrement supérieur à 2p" (ils ont p - la longueur de la mantisse avec la partie entière, c'est-à-dire p = 24 pour float). Plus loin dans le texte, ils montrent que la probabilité d'erreur avec une telle stratégie est proche de zéro, mais toujours positive, et cela est confirmé par des expériences.



Néanmoins, il existe encore des cas où la valeur de m doit être prise encore plus, et le pire des cas n'est pas connu à l'avance. Des estimations théoriques du pire des cas existent [1, section 12.7.2], mais elles donnent des millions de bits impensables, ce qui n'est pas bon. Voici un tableau des travaux cités (c'est pour la fonction exp (x) sur l'intervalle de -ln (2) à ln (2)):

Deuxième "béquille"

En pratique, m ne sera pas si terriblement grand. Et pour déterminer le pire des cas, une deuxième «béquille» est appliquée, qui est appelée «pré-calcul exhaustif». Pour le type de données float (32 bits), si la fonction f a un argument (x), alors nous pouvons facilement "exécuter" toutes les valeurs possibles de x. Le problème ne se posera qu'avec des fonctions qui ont plus d'un argument (parmi lesquelles pow (x, y)), pour lesquelles nous ne pourrions penser à rien de tel. Après avoir vérifié toutes les valeurs possibles de x, nous calculons notre constante m pour chaque fonction f (x) et pour chaque mode d'arrondi. Ensuite, les algorithmes de calcul qui doivent être implémentés dans le matériel sont conçus pour fournir une précision de 2 ^m . L'arrondi f (x) est alors garanti correct dans tous les cas.



Pour le type double (64 bits), une simple énumération est presque impossible. Cependant, ils trient! Mais comment? La réponse est donnée dans [4]. Je vais vous en parler très brièvement.



Le domaine de la fonction f (x) est divisé en très petits segments de sorte qu'à l'intérieur de chaque segment, il est possible de remplacer f (x) par une fonction linéaire de la forme b-ax (les coefficients a et b sont, bien entendu, différents pour des segments différents). La taille de ces segments est calculée analytiquement de sorte qu'une telle fonction linéaire serait en effet presque indiscernable de l'original dans chaque segment.



Ensuite, après quelques opérations de mise à l'échelle et de décalage, nous arrivons au problème suivant: une ligne droite b-ax peut-elle aller "assez près" d'un point entier?



Il s'avère qu'il est relativement facile de répondre par oui ou par non. Autrement dit, "oui" - si des points potentiellement dangereux sont proches d'une ligne droite, et "non" - si aucun point de ce type, en principe, ne peut se rapprocher de la ligne. La beauté de cette méthode est que la réponse «non» en pratique est obtenue dans l'écrasante majorité des cas, et la réponse «oui», rarement obtenue, vous oblige à parcourir le segment avec une force brute pour déterminer quels points spécifiques étaient critiques.



Cependant, l'itération sur les arguments de f (x) est réduite de nombreuses fois et permet de détecter des points de rupture pour des nombres comme double (binaire64) et long double (80 bits!). Cela se fait sur des supercalculateurs et, bien sûr, des cartes vidéo ... pendant leur temps libre de l'exploitation minière. Cependant, personne ne sait encore quoi faire avec le type de données binary128. Permettez-moi de vous rappeler que la partie fractionnaire de la mantisse de ces nombres est de 112 bits . Par conséquent, dans la littérature étrangère sur ce sujet jusqu'à présent, on ne peut trouver que des raisonnements semi-philosophiques commençant par «nous espérons ...» («nous espérons ...»).



Le détail de la méthode, qui permet de déterminer rapidement le passage d'une ligne à proximité de points entiers, est ici inapproprié. Pour ceux qui souhaitent apprendre le processus, je recommande de regarder plus attentivement le problème de trouver la distance entre une ligne droite et Z ² , par exemple, dans l'article [5]. Il décrit un algorithme amélioré qui, en cours de construction, ressemble au célèbre algorithme d'Euclid pour trouver le plus grand diviseur commun. Je donnerai la même image de [4] et [5], qui montre la transformation ultérieure du problème:







il existe de nombreux tableaux contenant les pires cas d'arrondi à des intervalles différents pour chaque fonction transcendantale. Ils se trouvent dans [1 section 12.8.4] et dans [3, section 10.5.3.2], ainsi que dans des articles séparés, par exemple, dans [6].



Je vais donner quelques exemples en prenant des lignes aléatoires à partir de ces tables. Je souligne que ce ne sont pas les pires cas pour tous les x, mais seulement pour quelques petits intervalles, voir la source si cela vous intéresse.







Comment lire le tableau? La valeur x est spécifiée en double notation hexadécimale à virgule flottante. Tout d'abord, comme prévu, il y a un premier, puis 52 bits de la partie fractionnaire de la mantisse et la lettre P. Cette lettre signifie "multiplier par deux à une puissance" suivi d'un degré. Par exemple, P-23 signifie que la mantisse spécifiée doit être multipliée par 2 ^-23 .



De plus, imaginez que la fonction f (x) soit calculée avec une précision infinie et que les 53 premiers bits en soient coupés (sans arrondir!). Ce sont ces 53 bits (dont un jusqu'à la virgule) qui sont indiqués dans la colonne f (x). Les bits suivants sont indiqués dans la dernière colonne. Le signe "degré" de la séquence de bits dans la dernière colonne signifie le nombre de répétitions de bits, c'est-à-dire, par exemple, 10 ⁵³1011 signifie que vient d'abord le bit égal à 1, puis 53 zéros et ensuite 1011. Ensuite, les points de suspension, ce qui signifie que nous n'avons en général pas besoin du reste des bits.



De plus, c'est une question de technologie - nous connaissons les pires cas pour chaque intervalle d'une fonction prise séparément et nous pouvons choisir pour cet intervalle une telle approximation afin qu'elle couvre le pire des cas avec sa précision. Avec seulement des années de supercalculateur, il est possible de créer des implémentations matérielles rapides et précises de fonctions élémentaires. Le problème est petit: il reste à apprendre au moins aux développeurs de compilateurs à utiliser ces tables.

Pourquoi est-ce nécessaire?

Excellente question! Après tout, j'ai dit à plusieurs reprises ci-dessus que près de 100% des programmeurs n'ont pas besoin de connaître une fonction élémentaire avec une précision au dernier bit correctement arrondi (souvent ils n'ont même pas besoin de la moitié des bits), pourquoi les scientifiques conduisent-ils des supercalculateurs et compilent des tableaux pour résoudre un problème «inutile»?



Premièrement, le défi est fondamental. Il est plutôt intéressant de ne pas obtenir un arrondi exact pour un arrondi précis, mais en principe pour comprendre comment ce problème intéressant pourrait être résolu, quels secrets des mathématiques computationnelles sa solution nous révélera-t-elle? Comment ces secrets pourraient-ils être utilisés dans d'autres tâches? Sciences fondamentales - elles sont comme ça, vous pouvez faire une sorte de "non-sens" pendant des décennies, puis cent ans plus tard, grâce à ce "non-sens", une percée scientifique se produit dans un autre domaine.



Deuxièmement, la question de la portabilité du code. Si une fonction peut se permettre de gérer les derniers bits du résultat de la manière qu'elle souhaite, cela signifie que sur différentes plates-formes et sur différents compilateurs, des résultats légèrement différents peuvent être obtenus (même si l'erreur spécifiée). Dans certains cas, ce n'est pas important, mais dans certains cas, cela peut être significatif, en particulier lorsque le programme a une erreur qui apparaît sur une plate-forme, mais qui n'apparaît pas sur une autre plate-forme précisément à cause des différents bits du résultat. Mais pourquoi est-ce que je vous décris le mal de tête bien connu associé aux différents comportements des programmes? Vous savez tout cela sans moi. Ce serait formidable d'avoir un système mathématique qui fonctionne exactement de la même manière sur toutes les plates-formes, quel que soit son niveau de compilation. C'est ce que vous devez faire correctement arrondissez le dernier morceau.

Liste des sources

[1] Jean-Michel Muller, «Fonctions élémentaires: algorithmes et implémentation», 2016



[2] William Kahan, « A Logarithm Too Clever by Half », 2004



[3] Jean-Michel Muller, «Handbook of floating-point arithmetic» , 2018



[4] Vincent Lefèvre, Jean-Michel Muller, «Vers des transcendantaux correctement arrondis», TRANSACTIONS IEEE SUR ORDINATEURS, VOL. 47, NO. 11, NOVEMBRE 1998, p. 1235-1243



[5] Vincent Lefèvre. «Nouveaux résultats sur la distance entre un segment et Z ² ». Application à l'arrondi exact. 17e Symposium IEEE sur l'arithmétique informatique - Arith'17, juin 2005, Cape Cod, MA,

États-Unis. 68 à 75



[6] Vincent Lefèvre, Jean-Michel Muller, «Pires cas d'arrondi correct des fonctions élémentaires en double précision», Rapport de recherche (INSTITUT NA TIONAL DE RECHERCHE EN INFORMA TIQUE ET EN AUTOMA TIQUE) n˚4044 - Novembre 2000 - 19 pages.