Résultat en virgule flottante différent avec l'optimisation activée - bogue du compilateur?

Le code ci-dessous fonctionne sur Visual Studio 2008 avec et sans optimisation. Mais cela ne fonctionne que sur g ++ sans optimisation (O0).

#include <cstdlib>
#include <iostream>
#include <cmath>

double round(double v, double digit)
{
    double pow = std::pow(10.0, digit);
    double t = v * pow;
    //std::cout << "t:" << t << std::endl;
    double r = std::floor(t + 0.5);
    //std::cout << "r:" << r << std::endl;
    return r / pow;
}

int main(int argc, char *argv[])
{
    std::cout << round(4.45, 1) << std::endl;
    std::cout << round(4.55, 1) << std::endl;
}

La sortie doit être:

4.5
4.6

Mais g ++ avec optimisation ( O1- O3) affichera:

4.5
4.5

Si j'ajoute le volatilemot - clé avant t, cela fonctionne, alors pourrait-il y avoir une sorte de bogue d'optimisation?

Test sur g ++ 4.1.2 et 4.4.4.

Voici le résultat sur ideone: http://ideone.com/Rz937

Et l'option que je teste sur g ++ est simple:

g++ -O2 round.cpp

Le résultat le plus intéressant, même si j'active l' /fp:fastoption sur Visual Studio 2008, le résultat est toujours correct.

Autre question:

Je me demandais, devrais-je toujours activer l' -ffloat-storeoption?

Parce que la version g ++ que j'ai testée est livrée avec CentOS / Red Hat Linux 5 et CentOS / Redhat 6 .

J'ai compilé beaucoup de mes programmes sous ces plates-formes et je crains que cela ne cause des bogues inattendus dans mes programmes. Il semble un peu difficile d'étudier tout mon code C ++ et les bibliothèques utilisées s'ils ont de tels problèmes. Toute suggestion?

Quelqu'un est-il intéressé par pourquoi même /fp:fastactivé, Visual Studio 2008 fonctionne toujours? Il semble que Visual Studio 2008 est plus fiable à ce problème que g ++?

Les processeurs Intel x86 utilisent une précision étendue de 80 bits en interne, alors qu'ils ont doublenormalement une largeur de 64 bits. Différents niveaux d'optimisation affectent la fréquence à laquelle les valeurs en virgule flottante du processeur sont enregistrées dans la mémoire et ainsi arrondies de la précision 80 bits à la précision 64 bits.

Utilisez l' -ffloat-storeoption gcc pour obtenir les mêmes résultats en virgule flottante avec différents niveaux d'optimisation.

Vous pouvez également utiliser le long doubletype, qui est normalement large de 80 bits sur gcc pour éviter d'arrondir une précision de 80 bits à 64 bits.

man gcc dit tout:

   -ffloat-store
       Do not store floating point variables in registers, and inhibit
       other options that might change whether a floating point value is
       taken from a register or memory.

       This option prevents undesirable excess precision on machines such
       as the 68000 where the floating registers (of the 68881) keep more
       precision than a "double" is supposed to have.  Similarly for the
       x86 architecture.  For most programs, the excess precision does
       only good, but a few programs rely on the precise definition of
       IEEE floating point.  Use -ffloat-store for such programs, after
       modifying them to store all pertinent intermediate computations
       into variables.

Dans les versions x86_64, les compilateurs utilisent des registres SSE pour floatet doublepar défaut, de sorte qu'aucune précision étendue ne soit utilisée et que ce problème ne se produise pas.

gccL'option du compilateur-mfpmath contrôle cela.

La sortie doit être: 4.5 4.6 C'est ce que serait la sortie si vous aviez une précision infinie, ou si vous travailliez avec un périphérique qui utilisait une représentation en virgule flottante basée sur des décimales plutôt que sur des bases binaires. Mais tu ne l'es pas. La plupart des ordinateurs utilisent la norme binaire à virgule flottante IEEE.

Comme Maxim Yegorushkin l'a déjà noté dans sa réponse, une partie du problème est que votre ordinateur utilise en interne une représentation en virgule flottante de 80 bits. Ceci n'est cependant qu'une partie du problème. La base du problème est que tout nombre de la forme n.nn5 n'a pas de représentation flottante binaire exacte. Ces cas de coin sont toujours des nombres inexacts.

Si vous voulez vraiment que votre arrondi puisse arrondir de manière fiable ces cas d'angle, vous avez besoin d'un algorithme d'arrondi qui tient compte du fait que n.n5, n.nn5 ou n.nnn5, etc. (mais pas n.5) est toujours inexact. Trouvez le cas d'angle qui détermine si une valeur d'entrée est arrondie vers le haut ou vers le bas et renvoie la valeur arrondie vers le haut ou arrondi vers le bas en fonction d'une comparaison avec ce cas d'angle. Et vous devez faire attention à ce qu'un compilateur optimisant ne placera pas ce cas d'angle trouvé dans un registre de précision étendue.

Voir Comment Excel arrondit avec succès les nombres flottants même s'ils sont imprécis? pour un tel algorithme.

Ou vous pouvez simplement vivre avec le fait que les boîtiers d'angle tournent parfois de manière erronée.

Différents compilateurs ont des paramètres d'optimisation différents. Certains de ces paramètres d'optimisation plus rapides ne maintiennent pas de règles strictes en virgule flottante selon IEEE 754 . Visual Studio dispose d' un cadre spécifique, /fp:strict, /fp:precise, /fp:fastoù /fp:fastviole la norme sur ce qui peut être fait. Vous constaterez peut-être que cet indicateur contrôle l'optimisation dans ces paramètres. Vous pouvez également trouver un paramètre similaire dans GCC qui change le comportement.

Si tel est le cas, la seule chose qui diffère entre les compilateurs est que GCC rechercherait par défaut le comportement en virgule flottante le plus rapide sur des optimisations plus élevées, alors que Visual Studio ne change pas le comportement en virgule flottante avec des niveaux d'optimisation plus élevés. Ainsi, ce n'est peut-être pas nécessairement un bogue réel, mais le comportement prévu d'une option que vous ne saviez pas que vous allumiez.

À ceux qui ne peuvent pas reproduire le bogue: ne pas décommenter les stmts de débogage commentés, ils affectent le résultat.

Cela implique que le problème est lié aux instructions de débogage. Et il semble qu'il y ait une erreur d'arrondi causée par le chargement des valeurs dans les registres pendant les instructions de sortie, c'est pourquoi d'autres ont trouvé que vous pouvez résoudre ce problème avec-ffloat-store

Autre question:

Je me demandais, devrais-je toujours activer l' -ffloat-storeoption?

Pour être irrévérencieux, il doit y avoir une raison que certains programmeurs ne tournent pas -ffloat-store, sinon l'option n'existerait pas ( de même, il doit y avoir une raison que certains programmeurs ne s'allument -ffloat-store). Je ne recommanderais pas de toujours l'activer ou de toujours l'éteindre. L'activer empêche certaines optimisations, mais sa désactivation permet le type de comportement que vous obtenez.

Mais, en général, il y a une certaine incohérence entre les nombres à virgule flottante binaire (comme l'ordinateur utilise) et les nombres décimaux à virgule flottante (avec lesquels les gens sont familiers), et cette discordance peut entraîner un comportement similaire à ce que vous obtenez (pour être clair, le comportement que vous obtenez n'est pas causé par cette incompatibilité, mais un comportement similaire peut être). Le fait est que, puisque vous avez déjà un peu d'imprécision en ce qui concerne la virgule flottante, je ne peux pas dire que -ffloat-storecela le rend meilleur ou pire.

Au lieu de cela, vous voudrez peut-être chercher d' autres solutions au problème que vous essayez de résoudre (malheureusement, Koenig ne pointe pas vers le papier réel, et je ne peux pas vraiment trouver un endroit "canonique" évident pour cela, donc je devra vous envoyer à Google ).

Si vous n'arrondissez pas à des fins de sortie, je regarderais probablement std::modf()(in cmath) et std::numeric_limits<double>::epsilon()(in limits). En réfléchissant à la round()fonction d' origine , je pense qu'il serait plus propre de remplacer l'appel à std::floor(d + .5)par un appel à cette fonction:

// this still has the same problems as the original rounding function
int round_up(double d)
{
    // return value will be coerced to int, and truncated as expected
    // you can then assign the int to a double, if desired
    return d + 0.5;
}

Je pense que cela suggère l'amélioration suivante:

// this won't work for negative d ...
// this may still round some numbers up when they should be rounded down
int round_up(double d)
{
    double floor;
    d = std::modf(d, &floor);
    return floor + (d + .5 + std::numeric_limits<double>::epsilon());
}

Une simple note: std::numeric_limits<T>::epsilon()est défini comme "le plus petit nombre ajouté à 1 qui crée un nombre différent de 1." Vous devez généralement utiliser un epsilon relatif (c'est-à-dire mettre à l'échelle epsilon pour tenir compte du fait que vous travaillez avec des nombres autres que "1"). La somme d, .5et std::numeric_limits<double>::epsilon()devrait être proche de 1, groupement ainsi que des moyens d'addition qui std::numeric_limits<double>::epsilon()seront de la bonne taille pour ce que nous faisons. Si quoi que ce soit, std::numeric_limits<double>::epsilon()sera trop grand (lorsque la somme des trois est inférieure à un) et peut nous amener à arrondir certains nombres alors que nous ne devrions pas.

De nos jours, vous devriez envisager std::nearbyint().

La réponse acceptée est correcte si vous compilez vers une cible x86 qui n'inclut pas SSE2. Tous les processeurs x86 modernes prennent en charge SSE2, donc si vous pouvez en profiter, vous devriez:

-mfpmath=sse -msse2 -ffp-contract=off

Décomposons cela.

-mfpmath=sse -msse2. Cela effectue l'arrondi à l'aide des registres SSE2, ce qui est beaucoup plus rapide que le stockage de chaque résultat intermédiaire en mémoire. Notez que c'est déjà la valeur par défaut sur GCC pour x86-64. Depuis le wiki GCC :

Sur les processeurs x86 plus modernes qui prennent en charge SSE2, la spécification des options du compilateur -mfpmath=sse -msse2garantit que toutes les opérations float et double sont effectuées dans les registres SSE et correctement arrondies. Ces options n'affectent pas l'ABI et doivent donc être utilisées autant que possible pour des résultats numériques prévisibles.

-ffp-contract=off. Cependant, contrôler l'arrondi ne suffit pas pour une correspondance exacte. Les instructions FMA (fused multiply-add) peuvent changer le comportement d'arrondi par rapport à ses homologues non fusionnés, nous devons donc le désactiver. Il s'agit de la valeur par défaut sur Clang, pas sur GCC. Comme expliqué par cette réponse :

Un FMA n'a qu'un seul arrondi (il conserve effectivement une précision infinie pour le résultat de multiplication temporaire interne), tandis qu'un ADD + MUL en a deux.

En désactivant FMA, nous obtenons des résultats qui correspondent exactement au débogage et à la publication, au prix de performances (et de précision). Nous pouvons toujours profiter des autres avantages de performance de SSE et AVX.

J'ai creusé davantage ce problème et je peux apporter plus de précisions. Premièrement, les représentations exactes de 4.45 et 4.55 selon gcc sur x84_64 sont les suivantes (avec libquadmath pour afficher la dernière précision):

float 32:   4.44999980926513671875
double 64:  4.45000000000000017763568394002504646778106689453125
doublex 80: 4.449999999999999999826527652402319290558807551860809326171875
quad 128:   4.45000000000000000000000000000000015407439555097886824447823540679418548304813185723105561919510364532470703125

float 32:   4.55000019073486328125
double 64:  4.54999999999999982236431605997495353221893310546875
doublex 80: 4.550000000000000000173472347597680709441192448139190673828125
quad 128:   4.54999999999999999999999999999999984592560444902113175552176459320581451695186814276894438080489635467529296875

Comme Maxim l'a dit ci-dessus, le problème est dû à la taille de 80 bits des registres FPU.

Mais pourquoi le problème ne se produit-il jamais sous Windows? sur IA-32, le FPU x87 était configuré pour utiliser une précision interne pour la mantisse de 53 bits (équivalent à une taille totale de 64 bits:) double. Pour Linux et Mac OS, la précision par défaut de 64 bits a été utilisée (équivalente à une taille totale de 80 bits :) long double. Le problème devrait donc être possible, ou pas, sur ces différentes plates-formes en changeant le mot de contrôle du FPU (en supposant que la séquence d'instructions déclencherait le bogue). Le problème a été signalé à gcc en tant que bogue 323 (lisez au moins le commentaire 92!).

Pour afficher la précision de la mantisse sous Windows, vous pouvez compiler ceci en 32 bits avec VC ++:

#include "stdafx.h"
#include <stdio.h>  
#include <float.h>  

int main(void)
{
    char t[] = { 64, 53, 24, -1 };
    unsigned int cw = _control87(0, 0);
    printf("mantissa is %d bits\n", t[(cw >> 16) & 3]);
}

et sur Linux / Cygwin:

#include <stdio.h>

int main(int argc, char **argv)
{
    char t[] = { 24, -1, 53, 64 };
    unsigned int cw = 0;
    __asm__ __volatile__ ("fnstcw %0" : "=m" (*&cw));
    printf("mantissa is %d bits\n", t[(cw >> 8) & 3]);
}

Notez qu'avec gcc, vous pouvez définir la précision FPU avec -mpc32/64/80, bien qu'elle soit ignorée dans Cygwin. Mais gardez à l'esprit que cela modifiera la taille de la mantisse, mais pas celle de l'exposant, laissant la porte ouverte à d'autres types de comportements différents.

Sur l'architecture x86_64, SSE est utilisé comme dit par tmandry , donc le problème ne se produira pas à moins que vous ne forciez l'ancien FPU x87 pour le calcul FP avec -mfpmath=387, ou à moins que -m32vous ne compiliez en mode 32 bits avec (vous aurez besoin du package multilib). Je pourrais reproduire le problème sous Linux avec différentes combinaisons de drapeaux et versions de gcc:

g++-5 -m32 floating.cpp -O1
g++-8 -mfpmath=387 floating.cpp -O1

J'ai essayé quelques combinaisons sur Windows ou Cygwin avec VC ++ / gcc / tcc mais le bogue n'est jamais apparu. Je suppose que la séquence d'instructions générée n'est pas la même.

Enfin, notez qu'une manière exotique d'éviter ce problème avec 4.45 ou 4.55 serait d'utiliser _Decimal32/64/128, mais le support est vraiment rare ... J'ai passé beaucoup de temps juste pour pouvoir faire un printf avec libdfp!

Personnellement, j'ai rencontré le même problème dans l'autre sens - de gcc à VS. Dans la plupart des cas, je pense qu'il vaut mieux éviter l'optimisation. Le seul moment où cela en vaut la peine, c'est lorsque vous avez affaire à des méthodes numériques impliquant de grands tableaux de données à virgule flottante. Même après le démontage, je suis souvent déçu par les choix des compilateurs. Très souvent, il est simplement plus facile d'utiliser les intrinsèques du compilateur ou simplement d'écrire l'assembly vous-même.