Pourquoi GCC utilise-t-il la multiplication par un nombre étrange pour implémenter la division d'un nombre entier ?

Question

J'ai lu des articles sur div y mul 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 ensuite générer du code en langage assembleur avec :

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

Mais en regardant ce qui est généré division.s il ne contient aucune opération de div ! À la place, il fait une sorte de magie noire avec des décalages 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

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

Answer 1

La division d'entiers est l'une des opérations arithmétiques les plus lentes que l'on puisse effectuer sur un processeur moderne, avec une latence pouvant atteindre des dizaines de cycles et un mauvais débit. (Pour x86, voir Tableaux d'instructions et guide de micro-architecture d'Agner Fog ).

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

Mise en œuvre du C / de cette manière plutôt qu'avec une séquence d'instructions multiples comprenant div est juste la façon par défaut de GCC de faire la division par des constantes. Elle ne nécessite pas d'optimisation entre les opérations et ne change rien, même pour le débogage. (En utilisant -Os pour une petite taille de code ne fait pas utiliser à GCC div mais ) Utiliser un inverse multiplicatif au lieu d'une division revient à utiliser lea au lieu de mul y add

Par conséquent, vous avez seulement tendance à voir div o idiv dans la sortie si 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 du code pour vous permettre de les générer vous-même (presque certainement inutile, à moins que vous ne travailliez avec un compilateur sans cervelle), voir libdivide .

Answer 2

Diviser par 5 revient à multiplier 1/5, ce qui revient à multiplier par 4/5 et à décaler de 2 bits vers la droite. La valeur concernée est CCCCCCCCCCCCCCCD en 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 quatre cinquièmes est 0.110011001100 récurrents - voir pourquoi ci-dessous). Je pense que vous pouvez prendre la suite ! Vous pouvez également consulter arithmétique à virgule fixe (mais notez qu'il est arrondi à un nombre entier à la fin).

Quant à savoir pourquoi, la multiplication est plus rapide que la division, et lorsque le diviseur est fixe, c'est une voie plus rapide.

Ver Multiplication réciproque, un tutoriel pour un article détaillé sur son fonctionnement, en expliquant en termes de points fixes. Il montre comment l'algorithme pour trouver la réciproque fonctionne, et comment gérer la division signée et le modulo.

Considérons un instant pourquoi 0.CCCCCCCC... (hex) ou 0.110011001100... Le binaire est 4/5. En divisant la représentation binaire par 4 (décalage de 2 places vers la droite), on obtient 0.001100110011... qui, par une inspection triviale, peut être ajouté à l'original pour obtenir 0.111111111111... qui est évidemment égal à 1, de la même manière que 0.9999999... en décimal est égal à un. Par conséquent, nous savons que x + x/4 = 1 donc 5x/4 = 1 , x=4/5 . Ceci est alors représenté par CCCCCCCCCCCCD en hex pour l'arrondi (car le chiffre binaire au-delà du dernier présent serait un 1 ).

Answer 3

En général, la multiplication est beaucoup plus rapide que la division. Donc, si on peut s'en sortir en multipliant par l'inverse, on peut accélérer considérablement la division par une constante.

Un problème est que nous ne pouvons pas représenter la réciproque exactement (à moins que la division ne soit une puissance de deux, mais dans ce cas, nous pouvons généralement convertir la division en un décalage de bits). Pour garantir des réponses correctes, nous devons donc veiller à ce que l'erreur dans notre réciproque n'entraîne 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 nombre entier de 64 bits par un nombre à virgule fixe de 0,64, nous obtenons un résultat de 64,64. Nous tronquons la valeur en un entier de 64 bits (en l'arrondissant effectivement à zéro), puis nous effectuons un nouveau décalage qui divise par quatre et tronque à nouveau. En regardant au niveau du bit, 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 faible pour ne pas pousser la réponse au-delà 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 de moins de 1/5 dans le résultat multiplié et décalé ne fera jamais passer le résultat au-delà d'une limite d'arrondi.

L'erreur dans notre constante est de (1/5) * 2 -64 . La valeur de i est inférieur à 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. -2 .

(1/5) * 2 -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 en arrondissant à partir de zéro, nous obtenons une erreur de (6/7) * 2. -64 . Après avoir multiplié par une valeur i d'un peu moins de 2 64 nous obtenons une erreur légèrement inférieure à 6/7 et, après avoir divisé par quatre, nous obtenons une erreur légèrement inférieure à 1,5/7, ce qui est supérieur à 1/7.

Ainsi, pour mettre en œuvre correctement la division par 7, nous devons multiplier par un nombre en virgule fixe de 0,65. Pour ce faire, il suffit de multiplier par les 64 bits inférieurs de notre nombre en virgule fixe, puis d'ajouter le nombre d'origine (ce qui peut déborder sur le bit de report) et enfin d'effectuer une rotation sur le report.

Answer 4

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 qui, 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 fixé à ceil(log2(d)), shpre est pre-shift (utilisé avant la multiplication) = e = nombre de bits zéro de queue dans d, shpost est post-shift (utilisé après la multiplication), prec est la précision = N - e = N - shpre. L'objectif est d'optimiser le calcul de n/d en utilisant un pré-décalage, une multiplication et un post-décalage.

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

La figure 4.2 et la figure 6.2 montrent comment le multiplicateur peut être réduit à un multiplicateur de 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 des processeurs modernes X86 et autres, le temps de multiplication est fixe, donc le pre-shift n'est pas utile 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 pre-shift pour les cibles X86.

Revenons à la figure 6.2. Le numérateur (dividende) pour mlow et mhigh peut être plus grand qu'un udword seulement quand le dénominateur (diviseur) > 2^(N-1) (quand == N => mlow = 2^(2N)), dans ce cas le remplacement optimisé pour 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 de N+1 bits, et deux divisions udword/uword peuvent être utilisées pour produire chaque valeur de N+1 bits (mlow ou mhigh). En utilisant 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 traité. Regardez comment j = i/7 est traité (ce qui devrait être le cas du multiplicateur à N+1 bits).

Sur la plupart des processeurs actuels, la multiplication a un timing fixe, donc un pre-shift 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 les 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
;       ...

Answer 5

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 méthode as-if règle.

• Le compilateur est autorisé à effectuer TOUTES les modifications tant qu'il ne change pas le comportement observable tel que 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 si beaucoup de programmeurs C le supposent). En général, il le fait parce que le compilateur veut optimiser les performances par rapport à l'approche directe (comme discuté en détail dans les autres réponses).
• Si, dans quelque circonstance que ce soit, le compilateur "optimise" un programme correct en quelque chose qui a un comportement observable différent, il s'agit d'un bogue de compilateur.
• Tout comportement non défini dans notre code (un dépassement d'entier signé est un exemple classique) et ce contrat est nul.