Division en virgule flottante

Dans mon code, il y a des endroits où je veux m’assurer qu’une division de 2 nombres à virgule flottante arbitraire (précision simple 32 bits) ne débordera pas. La cible / le compilateur ne garantit pas (assez explicitement) la manipulation agréable de -INF / INF et (ne garantit pas totalement IEEE 754 pour les valeurs exceptionnelles – (éventuellement non définie) – et la cible peut changer). De plus, je ne peux pas enregistrer d’hypothèses de sauvegarde sur les entrées pour ces quelques emplacements spéciaux et je suis lié aux bibliothèques standard C90.

J’ai lu ce que tous les informaticiens devraient savoir sur l’arithmétique en virgule flottante, mais pour être honnête, je suis un peu perdu.

Donc … Je veux demander à la communauté si le code suivant ferait l’affaire, et s’il existe des moyens meilleurs / plus rapides / correcteurs / correcteurs de le faire:

#define SIGN_F(val) ((val >= 0.0f)? 1.0f : -1.0f) float32_t safedivf(float32_t num, float32_t denum) { const float32_t abs_denum = fabs(denum); if((abs_denum < 1.0f) && ((abs_denum * FLT_MAX) <= (float32_t)fabs(num)) return SIGN_F(denum) * SIGN_F(num) * FLT_MAX; else return num / denum; } 

Edit: Modifié ((abs_denum * FLT_MAX) < (float32_t)fabs(num)) en ((abs_denum * FLT_MAX) <= (float32_t)fabs(num)) comme recommandé par Pascal Cuoq.

Dans ((abs_denum * FLT_MAX) < (float32_t)fabs(num) , le produit abs_denum * FLT_MAX peut être arrondi à la valeur fabs(num) . Cela ne signifie pas que num / denum n'est pas mathématiquement supérieur à FLT_MAX , et vous devriez craindre que cela cause le débordement que vous souhaitez éviter. Vous feriez mieux de remplacer ce < par <= .


Pour une solution alternative, si un type double est disponible et qu'il est plus large que float , il peut être plus économique de calculer (double)num/(double)denum . Si float est binary32ish et double est binary64ish, la division ne peut déborder que si denum est (a) zéro (et si denum est zéro, votre code est également problématique).

 double dbl_res = (double)num/(double)denum; float res = dbl_res < -FLT_MAX ? -FLT_MAX : dbl_res > FLT_MAX ? FLT_MAX : (float)dbl_res; 

Désolé, je ne peux pas commenter en raison de ma mauvaise réputation, mais vous pouvez essayer d’extraire les exposants et les mantisses de num et denum et vous assurer de cette condition:

 ((exp(num) - exp (denum)) > max_exp) && (mantissa(num) >= mantissa(denum)) 

Et en fonction du signe des entrées, générez le fichier INF correspondant.

Travaillez soigneusement avec num, denom lorsque le quotient est proche de FLT_MAX .

Les éléments suivants utilisent des tests inspirés de OP mais restnt à l’écart des résultats proches de FLT_MAX . Comme le souligne @Pascal Cuoq, le fait d’arrondir ne peut que pousser le résultat à la limite. Au lieu de cela, il utilise les seuils FLT_MAX/FLT_RADIX et FLT_MAX*FLT_RADIX .

En FLT_RADIX à l’échelle avec FLT_RADIX , généralement 2, le code doit toujours obtenir des résultats exacts. Arrondir dans n’importe quel mode d’arrondi ne devrait pas contaminer le résultat.

En termes de rapidité, le «bon chemin», c’est-à-dire lorsque les résultats ne sont certainement pas excessifs, devrait constituer un calcul rapide. Encore faut-il faire des tests unitaires, mais les commentaires devraient fournir l’essentiel de cette approche.

 static int SD_Sign(float x) { if (x > 0.0f) return 1; if (x < 0.0f) return -1; if (atan2f(x, -1.0f) > 0.0f) return 1; return -1; } static float SD_Overflow(float num, float denom) { return SD_Sign(num) * SD_Sign(denom) * FLT_MAX; } float safedivf(float num, float denom) { float abs_denom = fabsf(denom); // If |quotient| > |num| if (abs_denom < 1.0f) { float abs_num = fabsf(num); // If |num/denom| > FLT_MAX/2 --> quotient is very large or overflows // This computation is safe from rounding and overflow. if (abs_num > FLT_MAX / FLT_RADIX * abs_denom) { // If |num/denom| >= FLT_MAX*2 --> overflow // This also catches denom == 0.0 if (abs_num / FLT_RADIX >= FLT_MAX * abs_denom) { return SD_Overflow(num, denom); } // At this point, quotient must be in or near range FLT_MAX/2 to FLT_MAX*2 // Scale parameters so quotient is a FLT_RADIX * FLT_RADIX factor smaller. if (abs_num > 1.0) { abs_num /= FLT_RADIX * FLT_RADIX; } else { abs_denom *= FLT_RADIX * FLT_RADIX; } float quotient = abs_num / abs_denom; if (quotient > FLT_MAX / (FLT_RADIX * FLT_RADIX)) { return SD_Overflow(num, denom); } } } return num / denom; } 

Le SIGN_F() doit prendre en compte au +0.0 est +0.0 ou +0.0 . Diverses méthodes mentionnées par @Pascal Cuoq dans un commentaire:

  1. copysign() ou signbit()
  2. Utilisez un syndicat

De plus, certaines fonctions, lorsqu’elles sont bien implémentées, se différencient sur +/- zéro comme atan2f(zero, -1.0) et sprintf(buffer, "%+f", zero) .

Remarque: Utilisez float vs float32_t pour plus de simplicité.
Remarque: utilisez peut-être fabsf() plutôt que fabs() .
Mineur: Suggérez un dénominateur au lieu du denum .

Pour éviter les virages avec arrondissement et ce qui ne l’est pas, vous pouvez masser l’exposant sur le diviseur – avec frexp () et ldexp () – et vous demander si le résultat peut être réduit sans débordement. Ou frexp () les deux arguments, et effectue le travail de l’exposant à la main.