Implémentation de la division simple précision sous forme de multiplication double précision

Question

Pour un compilateur C99 mettant en œuvre l’arithmétique IEEE 754 exacte, existe-t-il des valeurs de f , divisor de type float telles que f / divisor != (float)(f * (1.0 / divisor)) ?

EDIT: Par «implémentation exacte de l’arithmétique IEEE 754», j’entends un compilateur qui définit à juste titre FLT_EVAL_METHOD sur 0.

Le contexte

Un compilateur CA qui fournit une virgule flottante conforme à la norme IEEE 754 ne peut remplacer une division simple précision par une constante par une multiplication simple par l’inverse que si ledit inverse est lui-même représentable exactement comme un float .

En pratique, cela ne se produit que pour des puissances de deux. Donc, un programmeur, Alex, peut être sûr que f / 2.0f sera compilé comme si c’était f * 0.5f , mais s’il est acceptable qu’Alex multiplie par 0.10f au lieu de diviser par 10, Alex devrait le dire par en écrivant la multiplication dans le programme, ou en utilisant une option de compilation, telle que -ffast-math de GCC.

Cette question concerne la transformation d’une division simple précision en une multiplication double précision. Le résultat est-il toujours arrondi? Y a-t-il une chance que cela soit moins cher, et donc une optimisation que les compilateurs pourraient faire (même sans -ffast-math )?

J’ai comparé (float)(f * 0.10) et f / 10.0f pour toutes les valeurs à simple précision de f comsockets entre 1 et 2, sans trouver aucun contre-exemple. Cela devrait couvrir toutes les divisions des float normaux produisant un résultat normal.

Ensuite, j’ai généralisé le test à tous les diviseurs avec le programme ci-dessous:

 #include  #include  #include  int main(void){ for (float divisor = 1.0; divisor != 2.0; divisor = nextafterf(divisor, 2.0)) { double factor = 1.0 / divisor; // double-precision inverse for (float f = 1.0; f != 2.0; f = nextafterf(f, 2.0)) { float cr = f / divisor; float opt = f * factor; // double-precision multiplication if (cr != opt) printf("For divisor=%a, f=%a, f/divisor=%a but (float)(f*factor)=%a\n", divisor, f, cr, opt); } } } 

L’espace de recherche est juste assez grand pour rendre cela intéressant (2 46 ). Le programme est en cours d’exécution. Quelqu’un peut-il me dire s’il imprimera quelque chose, peut-être avec une explication, pourquoi ou pourquoi pas, avant la fin?

Votre programme n’imprimera rien, en supposant que le mode d’arrondissement est égal à égal. L’essence de l’argument est la suivante:

Nous supposons que f et divisor sont compris entre 1.0 et 2.0 . Donc f = a / 2^23 et divisor = b / 2^23 pour certains entiers a et b compris dans l’intervalle [2^23, 2^24) . Le divisor = 1.0 cas divisor = 1.0 n’est pas intéressant, nous pouvons donc supposer que b > 2^23 .

Le seul moyen que (float)(f * (1.0 / divisor)) puisse donner un résultat erroné serait que la valeur exacte f / divisor soit aussi proche de la moitié du cas (c’est-à-dire un nombre situé exactement à mi-chemin entre deux floats) que les erreurs accumulées dans l’expression f * (1.0 / divisor) nous poussent de l’ autre côté de ce cas à mi-chemin par rapport à la valeur vraie.

Mais cela ne peut pas arriver. Par souci de simplicité, supposons d’abord que f >= divisor , de sorte que le quotient exact soit dans [1.0, 2.0) . Maintenant, tout cas à mi-parcours pour une précision simple dans l’intervalle [1.0, 2.0) a la forme c / 2^24 pour un entier impair c avec 2^24 < c < 2^25 . La valeur exacte de f / divisor est a / b , de sorte que la valeur absolue de la différence f / divisor - c / 2^24 est délimitée par 1 / (2^24 b) , donc au moins 1 / 2^48 (depuis b < 2^24 ). Nous sums donc à plus de 16 ulps en double précision dans tous les cas, et il devrait être facile de montrer que l'erreur dans le calcul en double précision ne peut jamais dépasser 16 ulps. (Je n'ai pas fait le calcul, mais je suppose qu'il est facile d'afficher une limite supérieure de 3 ulps sur l'erreur.)

Donc, le f / divisor ne peut pas être assez proche d'un cas à mi-chemin pour créer des problèmes. Notez que f / divisor ne peut pas être un cas à mi-chemin exact non plus: puisque c est impair, c et 2^24 sont relativement premiers, ainsi le seul moyen d’avoir c / 2^24 = a / b est si b est un multiple de 2^24 . Mais b est dans la plage (2^23, 2^24) , donc ce n'est pas possible.

Le cas où f < divisor est similaire: les cas intermédiaires ont alors la forme c / 2^25 et l’argument analogue montre que abs(f / divisor - c / 2^25) est supérieur à 1 / 2^49 , ce qui encore une fois. nous donne une marge de 16 ulps en double précision pour jouer.

Ce n’est certainement pas possible si des modes d’arrondi sans défaut sont possibles. Par exemple, en remplaçant 3.0f / 3.0f par 3.0f * C , une valeur de C inférieure à l’inverse exact donnerait un résultat erroné dans les modes d’arrondi vers le bas ou vers zéro, alors qu’une valeur de C supérieure à l’inverse exacte produirait donner un résultat erroné pour le mode d’arrondi à la hausse.

Il est moins clair pour moi si ce que vous recherchez est possible si vous vous limitez au mode d’arrondi par défaut. Je vais y réfléchir et réviser cette réponse si je trouve quelque chose.

Une recherche aléatoire a donné un exemple.

On dirait que lorsque le résultat est un nombre “dénormal / subnormal”, l’inégalité est possible. Mais alors, peut-être que ma plateforme n’est pas conforme à la norme IEEE 754?

 f 0x1.7cbff8p-25 divisor -0x1.839p+116 q -0x1.f8p-142 q2 -0x1.f6p-142 int MyIsFinite(float f) { union { float f; unsigned char uc[sizeof (float)]; unsigned long ul; } x; xf = f; return (x.ul & 0x7F800000L) != 0x7F800000L; } float floatRandom() { union { float f; unsigned char uc[sizeof (float)]; } x; do { size_t i; for (i=0; i 

Version PC Eclipse: Juno Service Release 2 Numéro de construction: 20130225-0426