¿Por qué GCC usa la multiplicación por un número extraño en la implementación de la división de enteros?

He estado leyendo sobre divy muloperaciones de montaje, y decidí ver en acción al escribir un programa sencillo en C:

File 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;
}

Y luego generar código de lenguaje ensamblador con:

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

Pero mirando el division.sarchivo generado , ¡no contiene ninguna operación div! En cambio, hace algún tipo de magia negra con pequeños cambios y números mágicos. Aquí hay un fragmento de código que calcula 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 está pasando aqui? ¿Por qué GCC no usa div en absoluto? ¿Cómo genera este número mágico y por qué todo funciona?

La división de enteros es una de las operaciones aritméticas más lentas que puede realizar en un procesador moderno, con una latencia de hasta docenas de ciclos y un rendimiento deficiente. (Para x86, consulte las tablas de instrucciones y la guía de microarquitectura de Agner Fog ).

Si conoce el divisor con anticipación, puede evitar la división reemplazándola con un conjunto de otras operaciones (multiplicaciones, sumas y cambios) que tengan el efecto equivalente. Incluso si se necesitan varias operaciones, a menudo sigue siendo mucho más rápido que la división entera en sí.

Implementar el /operador C de esta manera en lugar de con una secuencia de múltiples instrucciones involucradas dives solo la forma predeterminada de GCC de hacer la división por constantes. No requiere optimización en todas las operaciones y no cambia nada, incluso para la depuración. (Sin embargo, el uso -Ospara un tamaño de código pequeño hace que GCC lo use div). Usar un inverso multiplicativo en lugar de división es como usar en lealugar de mulyadd

Como resultado, solo tiende a ver divo idiven la salida si el divisor no se conoce en tiempo de compilación.

Para obtener información sobre cómo el compilador genera estas secuencias, así como el código que le permite generarlas por usted mismo (es casi innecesario a menos que esté trabajando con un compilador de braindead), vea libdivide .

Dividir por 5 es lo mismo que multiplicar 1/5, que es nuevamente lo mismo que multiplicar por 4/5 y desplazar a la derecha 2 bits. El valor en cuestión está CCCCCCCCCCCCCCCDen hexadecimal, que es la representación binaria de 4/5 si se coloca después de un punto hexadecimal (es decir, el binario para cuatro quintos es 0.110011001100recurrente; consulte a continuación para ver por qué). ¡Creo que puedes tomarlo desde aquí! Es posible que desee verificar la aritmética de punto fijo (aunque tenga en cuenta que se redondea a un número entero al final.

En cuanto a por qué, la multiplicación es más rápida que la división, y cuando el divisor es fijo, esta es una ruta más rápida.

Vea Multiplicación recíproca, un tutorial para una descripción detallada de cómo funciona, explicando en términos de punto fijo. Muestra cómo funciona el algoritmo para encontrar el recíproco y cómo manejar la división y el módulo con signo.

Consideremos por un minuto por qué 0.CCCCCCCC...(hexadecimal) o 0.110011001100...binario es 4/5. Divida la representación binaria entre 4 (cambie a la derecha 2 lugares), y obtendremos 0.001100110011...cuál por inspección trivial se puede agregar el original para obtener0.111111111111... , que obviamente es igual a 1, de la misma manera 0.9999999...en decimal es igual a uno. Por lo tanto, sabemos que x + x/4 = 1, por lo que 5x/4 = 1, x=4/5. Esto se representa como CCCCCCCCCCCCDen hexadecimal para redondear (ya que el dígito binario más allá del último presente sería a 1).

En general, la multiplicación es mucho más rápida que la división. Entonces, si podemos evitar la multiplicación por el recíproco, podemos acelerar significativamente la división por una constante

Una arruga es que no podemos representar el recíproco exactamente (a menos que la división fuera por una potencia de dos, pero en ese caso generalmente podemos convertir la división en un cambio de bits). Por lo tanto, para garantizar respuestas correctas, debemos tener cuidado de que el error en nuestro recíproco no cause errores en nuestro resultado final.

-3689348814741910323 es 0xCCCCCCCCCCCCCCCD, que es un valor de poco más de 4/5 expresado en 0.64 punto fijo.

Cuando multiplicamos un número entero de 64 bits por un número de punto fijo de 0.64 obtenemos un resultado de 64.64. Truncamos el valor a un entero de 64 bits (redondeándolo efectivamente a cero) y luego realizamos un cambio adicional que se divide entre cuatro y nuevamente se trunca Al observar el nivel de bits, está claro que podemos tratar ambas truncaciones como un solo truncamiento.

Esto claramente nos da al menos una aproximación de la división por 5, pero ¿nos da una respuesta exacta correctamente redondeada hacia cero?

Para obtener una respuesta exacta, el error debe ser lo suficientemente pequeño como para no empujar la respuesta sobre un límite de redondeo.

La respuesta exacta a una división por 5 siempre tendrá una parte fraccionaria de 0, 1/5, 2/5, 3/5 o 4/5. Por lo tanto, un error positivo de menos de 1/5 en el resultado multiplicado y desplazado nunca empujará el resultado sobre un límite de redondeo.

El error en nuestra constante es (1/5) * 2 ^-64 . El valor de i es menor que 2 ^64, por lo que el error después de multiplicar es menor que 1/5. Después de la división por 4, el error es menor que (1/5) * 2 ⁻² .

(1/5) * 2 ⁻² <1/5 por lo que la respuesta siempre será igual a hacer una división exacta y redondear hacia cero.

Lamentablemente, esto no funciona para todos los divisores.

Si intentamos representar 4/7 como un número de punto fijo de 0.64 con redondeo desde cero, terminamos con un error de (6/7) * 2 ^-64 . Después de multiplicar por un valor i de poco menos de 2 ⁶⁴ , terminamos con un error justo debajo de 6/7 y después de dividir entre cuatro terminamos con un error de poco menos de 1.5 / 7 que es mayor que 1/7.

Entonces, para implementar la división por 7 correctamente, debemos multiplicar por un número de punto fijo de 0.65. Podemos implementar eso multiplicando por los 64 bits más bajos de nuestro número de punto fijo, luego agregando el número original (esto puede desbordarse en el bit de acarreo) y luego haciendo una rotación a través del acarreo.

Aquí hay un enlace a un documento de un algoritmo que produce los valores y el código que veo con Visual Studio (en la mayoría de los casos) y que supongo que todavía se usa en GCC para la división de un entero variable por un entero constante.

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

En el artículo, una uword tiene N bits, una udword tiene 2N bits, n = numerador = dividendo, d = denominador = divisor, ℓ se establece inicialmente en ceil (log2 (d)), shpre es pre-shift (usado antes de multiplicar ) = e = número de bits cero finales en d, shpost es posterior al desplazamiento (utilizado después de la multiplicación), prec es precisión = N - e = N - shpre. El objetivo es optimizar el cálculo de n / d usando un pre-turno, multiplicación y post-turno.

Desplácese hasta la figura 6.2, que define cómo se genera un multiplicador de udwords (el tamaño máximo es N + 1 bits), pero no explica claramente el proceso. Explicaré esto a continuación.

La Figura 4.2 y la Figura 6.2 muestran cómo el multiplicador puede reducirse a un N bit o menos multiplicador para la mayoría de los divisores. La ecuación 4.5 explica cómo se obtuvo la fórmula utilizada para tratar con multiplicadores de N + 1 bit en las figuras 4.1 y 4.2.

En el caso de los procesadores X86 modernos y otros, el tiempo de multiplicación es fijo, por lo que el cambio previo no ayuda en estos procesadores, pero sí ayuda a reducir el multiplicador de N + 1 bits a N bits. No sé si GCC o Visual Studio han eliminado el cambio previo para los objetivos X86.

Volviendo a la Figura 6.2. El numerador (dividendo) para mlow y mhigh puede ser mayor que una udword solo cuando denominador (divisor)> 2 ^ (N-1) (cuando ℓ == N => mlow = 2 ^ (2N)), en este caso el El reemplazo optimizado para n / d es una comparación (si n> = d, q = 1, sino q = 0), por lo que no se genera un multiplicador. Los valores iniciales de mlow y mhigh serán N + 1 bits, y se pueden usar dos divisiones udword / uword para producir cada valor de N + 1 bit (mlow o mhigh). Usando X86 en modo de 64 bits como ejemplo:

; 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

Puedes probar esto con GCC. Ya has visto cómo se maneja j = i / 5. Eche un vistazo a cómo se maneja j = i / 7 (que debería ser el caso del multiplicador de N + 1 bit).

En la mayoría de los procesadores actuales, multiplicar tiene un tiempo fijo, por lo que no es necesario un cambio previo. Para X86, el resultado final es una secuencia de dos instrucciones para la mayoría de los divisores, y una secuencia de cinco instrucciones para divisores como 7 (para emular un multiplicador de N + 1 bit como se muestra en la ecuación 4.5 y la figura 4.2 del archivo pdf). Código de ejemplo 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
;       ...

Contestaré desde un ángulo ligeramente diferente: porque está permitido hacerlo.

C y C ++ se definen contra una máquina abstracta. El compilador transforma este programa en términos de la máquina abstracta en máquina concreta siguiendo la regla as-if .

• El compilador puede realizar CUALQUIER cambio siempre que no cambie el comportamiento observable según lo especificado por la máquina abstracta. No hay una expectativa razonable de que el compilador transformará su código de la manera más directa posible (incluso cuando muchos programadores de C asumen eso). Por lo general, hace esto porque el compilador quiere optimizar el rendimiento en comparación con el enfoque directo (como se explica en detalle en las otras respuestas).
• Si en alguna circunstancia el compilador "optimiza" un programa correcto para algo que tiene un comportamiento observable diferente, es un error del compilador.
• Cualquier comportamiento indefinido en nuestro código (desbordamiento de enteros firmados es un ejemplo clásico) y este contrato es nulo.