Pourquoi pow (int, int) est-il si lent?

J’ai travaillé sur quelques exercices de projet Euler pour améliorer mes connaissances en C ++.

J’ai écrit la fonction suivante:

int a = 0,b = 0,c = 0; for (a = 1; a <= SUMTOTAL; a++) { for (b = a+1; b <= SUMTOTAL-a; b++) { c = SUMTOTAL-(a+b); if (c == sqrt(pow(a,2)+pow(b,2)) && b < c) { std::cout << "a: " << a << " b: " << b << " c: "<< c << std::endl; std::cout << a * b * c << std::endl; } } } 

Cela se calcule en 17 millisecondes.

Cependant, si je change de ligne

 if (c == sqrt(pow(a,2)+pow(b,2)) && b < c) 

à

 if (c == sqrt((a*a)+(b*b)) && b < c) 

le calcul a lieu en 2 millisecondes. Y a-t-il des détails d’implémentation évidents de pow(int, int) qui me manquent qui rendent la première expression beaucoup plus lente?

    pow() fonctionne avec des nombres réels à virgule flottante et utilise la formule

     pow(x,y) = e^(y log(x)) 

    calculer x^y . Les int sont convertis pour double avant d’appeler pow . ( log est le logarithme naturel, basé sur e)

    x^2 utilisant pow() est donc plus lent que x*x .

    Modifier en fonction des commentaires pertinents

    • L’utilisation de pow même avec des exposants entiers peut donner des résultats incorrects ( PaulMcKenzie )
    • En plus d’utiliser une fonction mathématique avec double type, pow est un appel de fonction (alors que x*x n’est pas) ( jtbandes )
    • De nombreux compilateurs modernes vont en fait optimiser le pow avec des arguments entiers constants, mais il ne faut pas s’y fier.

    Vous avez choisi l’un des moyens les plus lents pour vérifier

     c*c == a*a + b*b // assuming c is non-negative 

    Cela se comstack en trois multiplications d’entiers (dont l’une peut être extraite de la boucle). Même sans pow() , vous convertissez toujours en double et prenez une racine carrée, ce qui est terrible pour le débit. (Et aussi la latence, mais la prédiction de twig + l’exécution spéculative sur les processeurs modernes signifie que la latence n’est pas un facteur ici).

    L’instruction SQRTSD d’Intel Haswell a un débit de 1 à 8-14 cycles ( source: tables d’instructions Agner Fog ), donc même si votre version de sqrt() maintient l’unité d’exécution FP sqrt saturée, elle rest environ 4 fois plus lente que gcc émettre (ci-dessous).


    Vous pouvez également optimiser la condition de la boucle pour sortir de la boucle lorsque la partie b < c de la condition devient fausse, de sorte que le compilateur doit uniquement effectuer une version de cette vérification.

     void foo_optimized() { for (int a = 1; a < = SUMTOTAL; a++) { for (int b = a+1; b < SUMTOTAL-ab; b++) { // int c = SUMTOTAL-(a+b); // gcc won't always transform signed-integer math, so this prevents hoisting (SUMTOTAL-a) :( int c = (SUMTOTAL-a) - b; // if (b >= c) break; // just changed the loop condition instead // the comstackr can hoist a*a out of the loop for us if (/* b < c && */ c*c == a*a + b*b) { // Just print a newline. std::endl also flushes, which bloats the asm std::cout << "a: " << a << " b: " << b << " c: "<< c << '\n'; std::cout << a * b * c << '\n'; } } } } 

    Cela comstack (avec gcc6.2 -O3 -mtune=haswell ) pour coder avec cette boucle interne. Voir le code complet sur l'explorateur du compilateur Godbolt .

     # a*a is hoisted out of the loop. It's in r15d .L6: add ebp, 1 # b++ sub ebx, 1 # c-- add r12d, r14d # ivtmp.36, ivtmp.43 # not sure what this is or why it's in the loop, would have to look again at the asm outside cmp ebp, ebx # b, _39 jg .L13 ## This is the loop-exit branch, not-taken until the end ## .L13 is the rest of the outer loop. ## It sets up for the next entry to this inner loop. .L8: mov eax, ebp # multiply a copy of the counters mov edx, ebx imul eax, ebp # b*b imul edx, ebx # c*c add eax, r15d # a*a + b*b cmp edx, eax # tmp137, tmp139 jne .L6 ## Fall-through into the cout print code when we find a match ## extremely rare, so should predict near-perfectly 

    Sur Intel Haswell, toutes ces instructions sont 1 uop chacune. (Et les paires cmp / jcc fusionnent dans des uops de comparaison et de twigment.) Donc, ce sont 10 uops de domaine fusionné, qui peuvent être émis à une itération par 2,5 cycles .

    Haswell exécute imul r32, r32 avec un débit d'une itération par horloge, de sorte que les deux multiplications à l'intérieur de la boucle interne ne saturent pas le port 1 à deux multiplications par 2,5c. Cela laisse de la place pour absorber les conflits de ressources inévitables entre ADD et SUB.

    Nous ne sums même pas proches des autres goulots d’étranglement des ports d’exécution, de sorte que le goulot d’étranglement frontal est le seul problème, qui devrait se produire à une itération par 2,5 cycles sur Intel Haswell et les versions ultérieures.

    Le déroulage en boucle pourrait aider ici à réduire le nombre de uops par contrôle. Par exemple, utilisez lea ecx, [rbx+1] pour calculer b + 1 pour la prochaine itération, afin que nous puissions imul ebx, ebx sans utiliser de MOV pour le rendre non destructif.


    Une réduction de la force est également possible : compte tenu de b*b nous pourrions essayer de calculer (b-1) * (b-1) sans IMUL. (b-1) * (b-1) = b*b - 2*b + 1 , donc peut-être que nous pouvons faire un lea ecx, [rbx*2 - 1] et ensuite soustraire cela de b*b . (Il n'y a pas de modes d'adressage qui soustraient au lieu d'append. Hmm, peut-être pourrions-nous garder -b dans un registre et compter jusqu'à zéro, donc nous pourrions utiliser lea ecx, [rcx + rbx*2 - 1] pour mettre à jour b*b dans ECX, donné -b dans EBX).

    À moins que vous ne gêniez réellement le débit IMUL, cela pourrait finir par prendre plus de poids et ne pas être une victoire. Il pourrait être amusant de voir à quel point un compilateur ferait bien avec cette réduction de la force dans le source C ++.


    Vous pouvez probablement également vectoriser ceci avec SSE ou AVX , en vérifiant 4 ou 8 valeurs b consécutives en parallèle. Comme les hits sont vraiment rares, il vous suffit de vérifier si l'un des 8 a eu un hit et de déterminer lequel était dans les rares cas où il y a eu une correspondance.

    Voir aussi le wiki des balises x86 pour plus d’optimisation.