Pourquoi GCC utilise-t-il la multiplication par un nombre étrange pour implémenter la division entière?

J'ai lu div et mulles opérations d'assemblage, et j'ai décidé de les voir en action en écrivant un programme simple en C:

Fichier division.c

#include <stdlib.h>
#include <stdio.h>

int main()
{
    size_t i = 9;
    size_t j = i / 5;
    printf("%zu\n",j);
    return 0;
}

Et puis générer du code en langage assembleur avec:

gcc -S division.c -O0 -masm=intel

Mais en regardant le division.sfichier généré , il ne contient aucune opération div! Au lieu de cela, il fait une sorte de magie noire avec un décalage de bits et des nombres magiques. Voici un extrait de code qui calcule i/5:

mov     rax, QWORD PTR [rbp-16]   ; Move i (=9) to RAX
movabs  rdx, -3689348814741910323 ; Move some magic number to RDX (?)
mul     rdx                       ; Multiply 9 by magic number
mov     rax, rdx                  ; Take only the upper 64 bits of the result
shr     rax, 2                    ; Shift these bits 2 places to the right (?)
mov     QWORD PTR [rbp-8], rax    ; Magically, RAX contains 9/5=1 now, 
                                  ; so we can assign it to j

Que se passe t-il ici? Pourquoi GCC n'utilise-t-il pas div du tout? Comment génère-t-il ce nombre magique et pourquoi tout fonctionne-t-il?

La division entière est l'une des opérations arithmétiques les plus lentes que vous puissiez effectuer sur un processeur moderne, avec une latence pouvant atteindre des dizaines de cycles et un mauvais débit. (Pour x86, voir les tableaux d'instructions et le guide microarch d'Agner Fog ).

Si vous connaissez le diviseur à l'avance, vous pouvez éviter la division en le remplaçant par un ensemble d'autres opérations (multiplications, ajouts et décalages) qui ont l'effet équivalent. Même si plusieurs opérations sont nécessaires, c'est souvent beaucoup plus rapide que la division entière elle-même.

L'implémentation de l' /opérateur C de cette façon au lieu d'une séquence multi-instructions impliquant divest juste la façon par défaut de GCC de faire la division par constantes. Il ne nécessite pas d'optimisation entre les opérations et ne change rien, même pour le débogage. (Utiliser -Ospour une petite taille de code permet cependant à GCC d'utiliser div.) L'utilisation d'un inverse multiplicatif au lieu de la division revient à utiliser à la leaplace de muletadd

Par conséquent, vous avez tendance à voir divou à afficher uniquement idivsi le diviseur n'est pas connu au moment de la compilation.

Pour plus d'informations sur la façon dont le compilateur génère ces séquences, ainsi que le code pour vous permettre de les générer vous-même (presque certainement inutile, sauf si vous travaillez avec un compilateur braindead ), voir libdivide .

Diviser par 5 équivaut à multiplier 1/5, ce qui équivaut à nouveau à multiplier par 4/5 et à décaler de 2 bits vers la droite. La valeur concernée est CCCCCCCCCCCCCCCDen hexadécimal, qui est la représentation binaire de 4/5 si elle est placée après un point hexadécimal (c'est-à-dire que le binaire pour les quatre cinquièmes est 0.110011001100récurrent - voir ci-dessous pour savoir pourquoi). Je pense que vous pouvez le prendre d'ici! Vous voudrez peut-être vérifier l' arithmétique à virgule fixe (bien que notez qu'il est arrondi à un entier à la fin.

Quant à savoir pourquoi, la multiplication est plus rapide que la division, et lorsque le diviseur est fixe, c'est un itinéraire plus rapide.

Voir Reciprocal Multiplication, un tutoriel pour une description détaillée de son fonctionnement, expliquant en termes de virgule fixe. Il montre comment fonctionne l'algorithme pour trouver la réciproque, et comment gérer la division signée et le modulo.

Voyons une minute pourquoi 0.CCCCCCCC...(hex) ou 0.110011001100...binaire est 4/5. Divisez la représentation binaire par 4 (décalage à droite de 2 places), et nous obtiendrons 0.001100110011...lequel, par inspection triviale, peut être ajouté l'original à obtenir 0.111111111111..., qui est évidemment égal à 1, de la même manière 0.9999999...en décimal est égal à un. Par conséquent, nous savons que x + x/4 = 1, si 5x/4 = 1, x=4/5. Ceci est alors représenté comme CCCCCCCCCCCCDen hexadécimal pour l'arrondi (comme le chiffre binaire au-delà du dernier présent serait a 1).

En général, la multiplication est beaucoup plus rapide que la division. Donc, si nous pouvons échapper à la multiplication par l'inverse, nous pouvons accélérer considérablement la division par une constante

Une ride est que nous ne pouvons pas représenter exactement l'inverse (à moins que la division soit par une puissance de deux, mais dans ce cas, nous pouvons généralement simplement convertir la division en un décalage de bits). Donc, pour garantir des réponses correctes, nous devons veiller à ce que l'erreur dans notre réciproque ne cause pas d'erreurs dans notre résultat final.

-3689348814741910323 est 0xCCCCCCCCCCCCCCCD qui est une valeur d'un peu plus de 4/5 exprimée en 0,64 point fixe.

Lorsque nous multiplions un entier 64 bits par un nombre à virgule fixe de 0,64, nous obtenons un résultat de 64,64. Nous tronquons la valeur à un entier 64 bits (l'arrondissant effectivement vers zéro), puis effectuons un autre décalage qui divise par quatre et tronque à nouveau.En regardant le niveau de bits, il est clair que nous pouvons traiter les deux troncatures comme une seule troncature.

Cela nous donne clairement au moins une approximation de la division par 5, mais cela nous donne-t-il une réponse exacte correctement arrondie vers zéro?

Pour obtenir une réponse exacte, l'erreur doit être suffisamment petite pour ne pas pousser la réponse au-dessus d'une limite d'arrondi.

La réponse exacte à une division par 5 aura toujours une partie fractionnaire de 0, 1/5, 2/5, 3/5 ou 4/5. Par conséquent, une erreur positive inférieure à 1/5 dans le résultat multiplié et décalé ne poussera jamais le résultat au-dessus d'une limite d'arrondi.

L'erreur dans notre constante est (1/5) * 2 ^-64 . La valeur de i est inférieure à 2 ^64, donc l'erreur après multiplication est inférieure à 1/5. Après la division par 4, l'erreur est inférieure à (1/5) * 2 ⁻² .

(1/5) * 2 ⁻² <1/5 donc la réponse sera toujours égale à faire une division exacte et à arrondir vers zéro.

Malheureusement, cela ne fonctionne pas pour tous les diviseurs.

Si nous essayons de représenter 4/7 comme un nombre à virgule fixe de 0,64 avec un arrondi à partir de zéro, nous nous retrouvons avec une erreur de (6/7) * 2 ^-64 . Après avoir multiplié par une valeur i d'un peu moins de 2 ^64, nous nous retrouvons avec une erreur un peu moins de 6/7 et après avoir divisé par quatre, nous nous retrouvons avec une erreur d'un peu moins de 1,5 / 7 qui est supérieure à 1/7.

Donc, pour implémenter correctement la division par 7, nous devons multiplier par un nombre à 0,65 point fixe. Nous pouvons implémenter cela en multipliant par les 64 bits inférieurs de notre nombre à virgule fixe, puis en ajoutant le nombre d'origine (cela peut déborder dans le bit de retenue) puis en effectuant une rotation par report.

Voici un lien vers un document d'un algorithme qui produit les valeurs et le code que je vois avec Visual Studio (dans la plupart des cas) et que je suppose est toujours utilisé dans GCC pour la division d'un entier variable par un entier constant.

http://gmplib.org/~tege/divcnst-pldi94.pdf

Dans l'article, un uword a N bits, un udword a 2N bits, n = numérateur = dividende, d = dénominateur = diviseur, ℓ est initialement défini sur ceil (log2 (d)), shpre est pré-décalage (utilisé avant multiplication ) = e = nombre de bits de fin à zéro dans d, shpost est post-shift (utilisé après multiplication), prec est précision = N - e = N - shpre. Le but est d'optimiser le calcul de n / d en utilisant un pré-décalage, une multiplication et un post-décalage.

Faites défiler jusqu'à la figure 6.2, qui définit comment un multiplicateur udword (la taille maximale est N + 1 bits), est généré, mais n'explique pas clairement le processus. Je vais l'expliquer ci-dessous.

Les figures 4.2 et 6.2 montrent comment le multiplicateur peut être réduit à un multiplicateur N bits ou moins pour la plupart des diviseurs. L'équation 4.5 explique comment la formule utilisée pour traiter les multiplicateurs N + 1 bits dans les figures 4.1 et 4.2 a été dérivée.

Dans le cas du X86 moderne et d'autres processeurs, le temps de multiplication est fixe, donc le pré-décalage n'aide pas sur ces processeurs, mais il aide toujours à réduire le multiplicateur de N + 1 bits à N bits. Je ne sais pas si GCC ou Visual Studio ont éliminé le pré-décalage pour les cibles X86.

Revenons à la figure 6.2. Le numérateur (dividende) pour mlow et mhigh ne peut être plus grand qu'un udword uniquement lorsque le dénominateur (diviseur)> 2 ^ (N-1) (lorsque ℓ == N => mlow = 2 ^ (2N)), dans ce cas, le le remplacement optimisé de n / d est une comparaison (si n> = d, q = 1, sinon q = 0), donc aucun multiplicateur n'est généré. Les valeurs initiales de mlow et mhigh seront N + 1 bits, et deux divisions udword / uword peuvent être utilisées pour produire chaque valeur N + 1 bits (mlow ou mhigh). Utiliser X86 en mode 64 bits comme exemple:

; upper 8 bytes of dividend = 2^(ℓ) = (upper part of 2^(N+ℓ))
; lower 8 bytes of dividend for mlow  = 0
; lower 8 bytes of dividend for mhigh = 2^(N+ℓ-prec) = 2^(ℓ+shpre) = 2^(ℓ+e)
dividend  dq    2 dup(?)        ;16 byte dividend
divisor   dq    1 dup(?)        ; 8 byte divisor

; ...
        mov     rcx,divisor
        mov     rdx,0
        mov     rax,dividend+8     ;upper 8 bytes of dividend
        div     rcx                ;after div, rax == 1
        mov     rax,dividend       ;lower 8 bytes of dividend
        div     rcx
        mov     rdx,1              ;rdx:rax = N+1 bit value = 65 bit value

Vous pouvez tester cela avec GCC. Vous avez déjà vu comment j = i / 5 est géré. Jetez un œil à la façon dont j = i / 7 est géré (qui devrait être le cas du multiplicateur N + 1 bits).

Sur la plupart des processeurs actuels, la multiplication a un timing fixe, donc un pré-décalage n'est pas nécessaire. Pour X86, le résultat final est une séquence de deux instructions pour la plupart des diviseurs, et une séquence de cinq instructions pour des diviseurs comme 7 (afin d'émuler un multiplicateur N + 1 bits comme indiqué dans l'équation 4.5 et la figure 4.2 du fichier pdf). Exemple de code X86-64:

;       rax = dividend, rbx = 64 bit (or less) multiplier, rcx = post shift count
;       two instruction sequence for most divisors:

        mul     rbx                     ;rdx = upper 64 bits of product
        shr     rdx,cl                  ;rdx = quotient
;
;       five instruction sequence for divisors like 7
;       to emulate 65 bit multiplier (rbx = lower 64 bits of multiplier)

        mul     rbx                     ;rdx = upper 64 bits of product
        sub     rbx,rdx                 ;rbx -= rdx
        shr     rbx,1                   ;rbx >>= 1
        add     rdx,rbx                 ;rdx = upper 64 bits of corrected product
        shr     rdx,cl                  ;rdx = quotient
;       ...

Je vais répondre sous un angle légèrement différent: Parce qu'il est permis de le faire.

C et C ++ sont définis par rapport à une machine abstraite. Le compilateur transforme ce programme en termes de machine abstraite en machine concrète en suivant la simulation règle .

• Le compilateur est autorisé à effectuer N'IMPORTE QUELLES modifications tant qu'il ne modifie pas le comportement observable spécifié par la machine abstraite. Il n'y a aucune attente raisonnable que le compilateur transforme votre code de la manière la plus simple possible (même lorsque beaucoup de programmeurs C le supposent). Habituellement, il le fait parce que le compilateur souhaite optimiser les performances par rapport à l'approche simple (comme discuté en détail dans les autres réponses).
• Si en aucun cas le compilateur "optimise" un programme correct en quelque chose qui a un comportement observable différent, c'est un bogue du compilateur.
• Tout comportement indéfini dans notre code (débordement d'entier signé est un exemple classique) et ce contrat est nul.