Pourquoi GCC n'optimise-t-il pas a * a * a * a * a * a à (a * a * a) * (a * a * a)?

Je fais une optimisation numérique sur une application scientifique. Une chose que j'ai remarquée est que GCC optimisera l'appel pow(a,2)en le compilant a*a, mais l'appel pow(a,6)n'est pas optimisé et appellera en fait la fonction de bibliothèque pow, ce qui ralentit considérablement les performances. (En revanche, Intel C ++ Compiler , exécutable icc, éliminera l'appel de bibliothèque pour pow(a,6).)

Ce que je suis curieux de savoir que lorsque je l' ai remplacé pow(a,6)avec l' a*a*a*a*a*aaide de GCC 4.5.1 et options « -O3 -lm -funroll-loops -msse4», il utilise 5 mulsdinstructions:

movapd  %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13

alors que si j'écris (a*a*a)*(a*a*a), ça produira

movapd  %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm14, %xmm13
mulsd   %xmm13, %xmm13

ce qui réduit le nombre d'instructions de multiplication à 3. icca un comportement similaire.

Pourquoi les compilateurs ne reconnaissent-ils pas cette astuce d'optimisation?

Parce que les mathématiques en virgule flottante ne sont pas associatives . La façon dont vous regroupez les opérandes en multiplication à virgule flottante a un effet sur la précision numérique de la réponse.

Par conséquent, la plupart des compilateurs sont très conservateurs quant à la réorganisation des calculs en virgule flottante, sauf s'ils peuvent être sûrs que la réponse restera la même ou si vous leur dites que vous ne vous souciez pas de la précision numérique. Par exemple: l' -fassociative-mathoption de gcc qui permet à gcc de réassocier des opérations en virgule flottante, ou même l' -ffast-mathoption qui permet des compromis encore plus agressifs entre précision et vitesse.

Lambdageek souligne correctement que, comme l'associativité n'est pas valable pour les nombres à virgule flottante, "l'optimisation" dea*a*a*a*a*ato(a*a*a)*(a*a*a)peut modifier la valeur. C'est pourquoi il est interdit par C99 (sauf si spécifiquement autorisé par l'utilisateur, via le drapeau du compilateur ou le pragma). En règle générale, l'hypothèse est que le programmeur a écrit ce qu'elle a fait pour une raison, et le compilateur doit respecter cela. Si tu veux(a*a*a)*(a*a*a), écris ça.

Cela peut cependant être pénible à écrire; pourquoi le compilateur ne peut-il pas simplement faire [ce que vous considérez être] la bonne chose lorsque vous l'utilisez pow(a,6)? Parce que ce serait la mauvaise chose à faire. Sur une plate-forme avec une bonne bibliothèque mathématique, pow(a,6)est nettement plus précis que l'un a*a*a*a*a*aou l' autre (a*a*a)*(a*a*a). Juste pour fournir quelques données, j'ai effectué une petite expérience sur mon Mac Pro, mesurant la pire erreur dans l'évaluation d'un ^ 6 pour tous les nombres flottants simple précision entre [1,2):

worst relative error using    powf(a, 6.f): 5.96e-08
worst relative error using (a*a*a)*(a*a*a): 2.94e-07
worst relative error using     a*a*a*a*a*a: 2.58e-07

L'utilisation powau lieu d'un arbre de multiplication réduit l'erreur liée par un facteur de 4 . Les compilateurs ne doivent pas (et ne le font généralement pas) faire des "optimisations" qui augmentent les erreurs à moins que l'utilisateur ne le permette (par exemple via -ffast-math).

Notez que GCC fournit __builtin_powi(x,n)une alternative à pow( ), qui devrait générer un arbre de multiplication en ligne. Utilisez-le si vous souhaitez échanger la précision contre les performances, mais ne souhaitez pas activer les mathématiques rapides.

Un autre cas similaire: la plupart des compilateurs n'optimiseront a + b + c + dpas (a + b) + (c + d)(c'est une optimisation car la deuxième expression peut être mieux canalisée) et l'évalueront comme donnée (c'est-à-dire comme (((a + b) + c) + d)). C'est aussi à cause des cas d'angle:

float a = 1e35, b = 1e-5, c = -1e35, d = 1e-5;
printf("%e %e\n", a + b + c + d, (a + b) + (c + d));

Cette sorties 1.000000e-05 0.000000e+00

Fortran (conçu pour le calcul scientifique) a un opérateur de puissance intégré, et pour autant que je sache, les compilateurs Fortran optimisent généralement l'augmentation des puissances entières d'une manière similaire à ce que vous décrivez. C / C ++ n'a malheureusement pas d'opérateur de puissance, seulement la fonction de bibliothèque pow(). Cela n'empêche pas les compilateurs intelligents de les traiter powspécialement et de les calculer plus rapidement pour des cas particuliers, mais il semble qu'ils le fassent moins souvent ...

Il y a quelques années, j'essayais de rendre plus pratique le calcul optimal des puissances entières, et j'ai proposé ce qui suit. C'est du C ++, pas du C cependant, et cela dépend toujours du fait que le compilateur soit quelque peu intelligent sur la façon d'optimiser / incorporer les choses. Quoi qu'il en soit, j'espère que cela vous sera utile dans la pratique:

template<unsigned N> struct power_impl;

template<unsigned N> struct power_impl {
    template<typename T>
    static T calc(const T &x) {
        if (N%2 == 0)
            return power_impl<N/2>::calc(x*x);
        else if (N%3 == 0)
            return power_impl<N/3>::calc(x*x*x);
        return power_impl<N-1>::calc(x)*x;
    }
};

template<> struct power_impl<0> {
    template<typename T>
    static T calc(const T &) { return 1; }
};

template<unsigned N, typename T>
inline T power(const T &x) {
    return power_impl<N>::calc(x);
}

_{Clarification pour les curieux: cela ne trouve pas le moyen optimal de calculer les puissances, mais comme la recherche de la solution optimale est un problème NP-complet et que cela ne vaut de toute façon que pour les petites puissances (par opposition à l'utilisation pow), il n'y a aucune raison de s'inquiéter avec le détail.}

Ensuite, utilisez-le simplement comme power<6>(a).

Cela facilite la saisie des pouvoirs (pas besoin d'énoncer 6 as avec des parens), et vous permet d'avoir ce type d'optimisation sans -ffast-mathau cas où vous auriez quelque chose qui dépend de la précision comme la sommation compensée (un exemple où l'ordre des opérations est essentiel) .

Vous pouvez aussi probablement oublier qu'il s'agit de C ++ et simplement l'utiliser dans le programme C (s'il compile avec un compilateur C ++).

J'espère que cela peut être utile.

ÉDITER:

Voici ce que j'obtiens de mon compilateur:

Pour a*a*a*a*a*a,

    movapd  %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0

Pour (a*a*a)*(a*a*a),

    movapd  %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm1, %xmm0
    mulsd   %xmm0, %xmm0

Pour power<6>(a),

    mulsd   %xmm0, %xmm0
    movapd  %xmm0, %xmm1
    mulsd   %xmm0, %xmm1
    mulsd   %xmm0, %xmm1

GCC n'optimise réellement a*a*a*a*a*aà (a*a*a)*(a*a*a)quand a est un entier. J'ai essayé avec cette commande:

$ echo 'int f(int x) { return x*x*x*x*x*x; }' | gcc -o - -O2 -S -masm=intel -x c -

Il y a beaucoup de drapeaux gcc mais rien d'extraordinaire. Ils signifient: Lire depuis stdin; utiliser le niveau d'optimisation d'O2; liste de langage d'assemblage de sortie au lieu d'un binaire; la liste doit utiliser la syntaxe du langage d'assemblage Intel; l'entrée est en langage C (généralement la langue est déduite de l'extension du fichier d'entrée, mais il n'y a pas d'extension de fichier lors de la lecture depuis stdin); et écrivez à stdout.

Voici la partie importante de la sortie. Je l'ai annoté avec quelques commentaires indiquant ce qui se passe dans le langage d'assemblage:

; x is in edi to begin with.  eax will be used as a temporary register.
mov  eax, edi  ; temp = x
imul eax, edi  ; temp = x * temp
imul eax, edi  ; temp = x * temp
imul eax, eax  ; temp = temp * temp

J'utilise le système GCC sur Linux Mint 16 Petra, un dérivé d'Ubuntu. Voici la version gcc:

$ gcc --version
gcc (Ubuntu/Linaro 4.8.1-10ubuntu9) 4.8.1

Comme d'autres affiches l'ont noté, cette option n'est pas possible en virgule flottante, car l'arithmétique en virgule flottante n'est pas associative.

Parce qu'un nombre à virgule flottante 32 bits - tel que 1,024 - n'est pas 1,024. Dans un ordinateur, 1.024 est un intervalle: de (1.024-e) à (1.024 + e), où "e" représente une erreur. Certaines personnes ne réalisent pas cela et croient également que * dans a * a signifie la multiplication de nombres de précision arbitraire sans qu'il y ait des erreurs attachées à ces nombres. La raison pour laquelle certaines personnes ne réalisent pas cela est peut-être les calculs mathématiques qu'elles ont exercés dans les écoles élémentaires: travailler uniquement avec des nombres idéaux sans erreurs attachées, et croire qu'il est OK d'ignorer simplement "e" tout en effectuant la multiplication. Ils ne voient pas le "e" implicite dans "float a = 1.2", "a * a * a" et les codes C similaires.

Si la majorité des programmeurs reconnaissent (et sont capables d'exécuter) l'idée que l'expression C a * a * a * a * a * a ne fonctionne pas réellement avec des nombres idéaux, le compilateur GCC serait alors GRATUIT pour optimiser "a * a * a * a * a * a "en dire" t = (a * a); t * t * t "qui nécessite un plus petit nombre de multiplications. Mais malheureusement, le compilateur GCC ne sait pas si le programmeur qui écrit le code pense que "a" est un nombre avec ou sans erreur. Et donc GCC ne fera que ce à quoi ressemble le code source - parce que c'est ce que GCC voit à l'œil nu.

... une fois que vous savez quel type de programmeur vous êtes, vous pouvez utiliser le commutateur "-ffast-math" pour dire à GCC que "Hé, GCC, je sais ce que je fais!". Cela permettra à GCC de convertir a * a * a * a * a * a en un autre morceau de texte - il semble différent de a * a * a * a * a * a - mais calcule toujours un nombre dans l'intervalle d'erreur de a * a * a * a * a * a. C'est OK, car vous savez déjà que vous travaillez avec des intervalles, pas des nombres idéaux.

Aucune affiche n'a encore mentionné la contraction des expressions flottantes (norme ISO C, 6.5p8 et 7.12.2). Si le FP_CONTRACTpragma est défini sur ON, le compilateur est autorisé à considérer une expression telle qu'une a*a*a*a*a*aopération unique, comme si elle était évaluée exactement avec un seul arrondi. Par exemple, un compilateur peut le remplacer par une fonction d'alimentation interne à la fois plus rapide et plus précise. Ceci est particulièrement intéressant car le comportement est en partie contrôlé par le programmeur directement dans le code source, tandis que les options du compilateur fournies par l'utilisateur final peuvent parfois être utilisées de manière incorrecte.

L'état par défaut du FP_CONTRACTpragma est défini par l'implémentation, de sorte qu'un compilateur est autorisé à effectuer de telles optimisations par défaut. Ainsi, le code portable qui doit suivre strictement les règles IEEE 754 doit le définir explicitement OFF.

Si un compilateur ne prend pas en charge ce pragma, il doit être prudent en évitant une telle optimisation, au cas où le développeur aurait choisi de le définir OFF.

GCC ne prend pas en charge ce pragma, mais avec les options par défaut, il le suppose ON; ainsi pour les cibles avec un FMA matériel, si l'on veut empêcher la transformation a*b+cen fma (a, b, c), il faut fournir une option telle que -ffp-contract=off(pour définir explicitement le pragma sur OFF) ou -std=c99(pour dire à GCC de se conformer à certains Version standard C, ici C99, suivez donc le paragraphe ci-dessus). Dans le passé, cette dernière option n'empêchait pas la transformation, ce qui signifie que GCC n'était pas conforme sur ce point: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=37845

Comme Lambdageek l'a souligné, la multiplication par flotteur n'est pas associative et vous pouvez obtenir moins de précision, mais aussi lorsque vous obtenez une meilleure précision, vous pouvez vous opposer à l'optimisation, car vous voulez une application déterministe. Par exemple, dans le client / serveur de simulation de jeu, où chaque client doit simuler le même monde que vous voulez que les calculs en virgule flottante soient déterministes.

Les fonctions de bibliothèque comme "pow" sont généralement soigneusement conçues pour produire le minimum d'erreur possible (dans le cas générique). Ceci est généralement réalisé en rapprochant les fonctions avec des splines (selon le commentaire de Pascal, l'implémentation la plus courante semble utiliser l' algorithme Remez )

fondamentalement l'opération suivante:

pow(x,y);

a une erreur inhérente d'environ la même ampleur que l'erreur dans une multiplication ou une division unique .

Alors que l'opération suivante:

float a=someValue;
float b=a*a*a*a*a*a;

a une erreur inhérente supérieure à 5 fois l'erreur d'une seule multiplication ou division (car vous combinez 5 multiplications).

Le compilateur doit faire très attention au type d'optimisation qu'il fait:

1. si l' optimisation pow(a,6)de a*a*a*a*a*ace peut améliorer les performances, mais aussi de réduire considérablement la précision des nombres à virgule flottante.
2. si l' optimisation a*a*a*a*a*a de pow(a,6)cela peut réellement réduire la précision parce que « a » est une valeur spéciale qui permet la multiplication sans erreur (une puissance de 2 ou d' un petit nombre entier)
3. si l' optimisation pow(a,6)à (a*a*a)*(a*a*a)ou (a*a)*(a*a)*(a*a)il peut encore être une perte de précision par rapport à la powfonction.

En général, vous savez que pour les valeurs arbitraires en virgule flottante, "pow" a une meilleure précision que n'importe quelle fonction que vous pourriez éventuellement écrire, mais dans certains cas spéciaux, les multiplications multiples peuvent avoir une meilleure précision et de meilleures performances, c'est au développeur de choisir ce qui est le plus approprié, commentant éventuellement le code afin que personne d'autre "optimise" ce code.

La seule chose qui ait du sens (opinion personnelle et apparemment un choix dans GCC sans optimisation particulière ou indicateur de compilateur) à optimiser devrait être de remplacer "pow (a, 2)" par "a * a". Ce serait la seule chose sensée qu'un vendeur de compilateur devrait faire.

Je ne m'attendais pas du tout à ce que ce cas soit optimisé. Cela ne peut pas être très souvent lorsqu'une expression contient des sous-expressions qui peuvent être regroupées pour supprimer des opérations entières. Je m'attendrais à ce que les rédacteurs de compilateurs investissent leur temps dans des domaines qui seraient plus susceptibles d'entraîner des améliorations notables, plutôt que de couvrir un cas de bord rarement rencontré.

J'ai été surpris d'apprendre des autres réponses que cette expression pouvait en effet être optimisée avec les bons commutateurs du compilateur. Soit l'optimisation est triviale, soit il s'agit d'un cas extrême d'une optimisation beaucoup plus courante, soit les rédacteurs du compilateur étaient extrêmement approfondis.

Il n'y a rien de mal à fournir des astuces au compilateur comme vous l'avez fait ici. C'est une partie normale et attendue du processus de micro-optimisation de réorganiser les instructions et les expressions pour voir quelles différences elles apporteront.

Bien que le compilateur puisse être justifié de considérer les deux expressions pour fournir des résultats incohérents (sans les commutateurs appropriés), il n'est pas nécessaire que vous soyez lié par cette restriction. La différence sera incroyablement minime - à tel point que si la différence vous tient à cœur, vous ne devriez pas utiliser l'arithmétique à virgule flottante standard en premier lieu.

Il y a déjà quelques bonnes réponses à cette question, mais dans un souci d'exhaustivité, je voulais souligner que la section applicable de la norme C est 5.1.2.2.3 / 15 (qui est la même que la section 1.9 / 9 du Norme C ++ 11). Cette section indique que les opérateurs ne peuvent être regroupés que s'ils sont réellement associatifs ou commutatifs.

gcc peut réellement faire cette optimisation, même pour les nombres à virgule flottante. Par exemple,

double foo(double a) {
  return a*a*a*a*a*a;
}

devient

foo(double):
    mulsd   %xmm0, %xmm0
    movapd  %xmm0, %xmm1
    mulsd   %xmm0, %xmm1
    mulsd   %xmm1, %xmm0
    ret

avec -O -funsafe-math-optimizations. Cette réorganisation viole IEEE-754, cependant, elle nécessite donc le drapeau.

Les entiers signés, comme Peter Cordes l'a souligné dans un commentaire, peuvent faire cette optimisation sans -funsafe-math-optimizationscar il tient exactement quand il n'y a pas de débordement et s'il y a débordement, vous obtenez un comportement indéfini. Vous obtenez donc

foo(long):
    movq    %rdi, %rax
    imulq   %rdi, %rax
    imulq   %rdi, %rax
    imulq   %rax, %rax
    ret

avec juste -O. Pour les entiers non signés, c'est encore plus facile car ils fonctionnent avec des puissances de mod de 2 et peuvent donc être réorganisés librement même en cas de débordement.