Est-ce qu'un réciproque à virgule flottante toujours aller-retour?

Question

Pour l'arithmétique IEEE-754, y a-t-il une garantie de 0 ou 1 unité à la dernière place pour la précision des réciproques? À partir de là, y a-t-il une garantie d'erreur sur l'inverse d'un réciproque?

Answer 1

[Tout ci-dessous suppose un format binaire IEEE 754 fixe, avec une forme de ronde à la plus proche que l'arrondissement de mode.]

Depuis réciproque (calculé comme 1/x) est une opération arithmétique de base, 1 est exactement représentable, et les opérations arithmétiques sont garantis correctement arrondis par le standard, le résultat réciproque est garanti à 0.5 unités à la dernière place de la vraie valeur. (Ceci est valable pour tout des opérations arithmétiques de base spécifiées dans la norme.)

L'inverse de la réciproque d'une valeur x ne sont pas garanties pour être égale à x, en général. Un exemple rapide avec le format IEEE 754 binary64:

>>> x = 1.8 
>>> 1.0/(1.0/x) 
1.7999999999999998 
>>> x == 1.0/(1.0/x) 
False 

Cependant, en supposant que le débordement et le sousécoulement sont évités, et que x est fini, non nulle, et sont vraies exactement représentable, les résultats suivants:

1. La valeur de 1.0/(1.0/x) sera différente de x de pas plus de 1 unité à la dernière place.

2. Si le significand de x (normalisé pour se situer dans la plage [1.0, 2.0) comme d'habitude) est plus petit que sqrt(2), l'réciproque-t aller-retour: 1.0/(1.0/x) == x.

Esquisse de la preuve: sans perte de généralité, on peut supposer que x est positif, et à l'échelle x par une puissance de deux pour qu'il se trouve dans la gamme [1.0, 2.0). Les résultats ci-dessus sont manifestement vrais dans le cas où x est une puissance exacte de deux, alors supposons que ce n'est pas le cas (cela sera utile dans la deuxième étape ci-dessous). La preuve ci-dessous est donnée pour la précision spécifique au format IEEE 754 binary64, mais elle s'adapte directement à tout format binaire IEEE 754.

Ecrire 1/x pour la vraie valeur de l'inverse, avant l'arrondissement, et laisser y être le (unique, il se trouve) le plus proche flotteur représentable binary64 à 1/x. Puis:

• depuis y est le plus proche flotteur 1/x, et les deux y et 1/x sont dans le binade [0.5, 1.0], où l'espacement entre les flotteurs successifs est exactement 2^-53, nous avons |y - 1/x| <= 2^-54. En fait, nous pouvons faire un peu mieux:

• nous avons en fait une inégalité stricte ci-dessus: |y - 1/x| < 2^-54. Si |y - 1/x| étaient exactement égal à 2^-54, alors 1/x serait exactement représentable en point flottant binaire de précision arbitraire (puisque les deux y et 2^-54 sont).Mais les seuls flottants binaires x pour lesquels 1/x est exactement représentable à une certaine précision sont des puissances de deux, et nous avons déjà exclu ce cas.

• si x < sqrt(2) puis 1/x > x/2, d'où (arrondi à la fois au flotteur représentable le plus proche), nous avons y >= x/2, donc x/y <= 2.

• maintenant x - 1/y = (y - 1/x) x/y et des limites sur |y - 1/x| et x/y (en supposant encore que x < sqrt(2)) nous obtenons |x - 1/y| < 2^-53. Il s'ensuit que x est le flottant représentable le plus proche de 1/y, 1/y arrondit à x et l'aller-retour réussit. Cela complète la preuve de la partie 2.

• dans le cas général x < 2, nous avons x/y < 4, donc |x - 1/y| < 2^-52. Cela fait 1/y au plus 1 ULP loin de x, qui complète la preuve d'une partie 1.

Voici une démonstration du seuil sqrt(2): en utilisant Python, nous prenons un million de flotteurs au hasard dans la gamme [1.0, 2.0), et identifier ceux qui ne font pas le tour de la réciproque. Tous les échantillons plus petits que sqrt(2) passent le roundtrip.

>>> import random 
>>> samples = [random.uniform(1.0, 2.0) for _ in range(10**6)] 
>>> bad = [x for x in samples if 1.0/(1.0/x) != x] 
>>> len(bad) 
171279 
>>> min(bad) 
1.4150519879892107 
>>> import math 
>>> math.sqrt(2) 
1.4142135623730951 

et une démonstration que l'erreur maximale est de pas plus de 1 ULP, en général (pour le format de binary64, dans le binade [1.0, 2.0), une unité à la dernière place est 2^-52) :

>>> samples = [random.uniform(1.0, 2.0) for _ in range(10**6)] 
>>> max(abs(x - 1.0/(1.0/x)) for x in samples) 
2.220446049250313e-16 
>>> 2**-52 
2.220446049250313e-16 

Voici un exemple avec le format IEEE 754 binary64 montre que la restriction à éviter soupassement est nécessaire:

>>> x = 1.3e308 
>>> x_roundtrip = 1.0/(1.0/x) 
>>> x.hex() 
'0x1.72409614c1e6ap+1023' 
>>> x_roundtrip.hex() 
'0x1.72409614c1e6cp+1023' 

Voici le résultat x_roundtrip diffère du b d'origine y deux unités à la dernière place, car 1/x était plus petit que le plus petit flottant représentable normal, et n'était donc pas représenté avec la même précision que x. Note finale: étant donné que la norme IEEE 754-2008 couvre également les types à virgule flottante décimale, je dois mentionner que la preuve ci-dessus est presque textuelle dans le cas décimal, établissant que pour les flottants de moins significatifs que sqrt(10), alors que pour les flottants décimaux généraux (évitant à nouveau le débordement et le débordement), nous ne sommes jamais éloignés de plus d'une unité à la dernière place. Cependant, il y a une certaine finesse théorique nécessaire pour montrer que l'inégalité clé |x - 1/y| < 1/2 10^(1-p) est toujours stricte: on finit par devoir montrer que la quantité 1 + 16 10^(2p-1) n'est jamais un nombre carré (ce qui est vrai, mais c'est sans doute hors de la portée de ce site inclure la preuve ici).

>>> from decimal import Decimal, getcontext 
>>> import random 
>>> getcontext().prec = 6 
>>> samples = [+Decimal(random.uniform(1.0, 10.0)) for _ in range(10**6)] 
>>> bad = [x for x in samples if 1/(1/x) != x] 
>>> min(bad) 
Decimal('3.16782')