¿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 acerca de las operaciones de ensamblaje div y mul , y decidí verlos en acción escribiendo un progtwig simple en C:

File division.c

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

Y luego generar el código de lenguaje ensamblador con:

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

¡Pero mirando el archivo generado de division.s , no contiene ninguna operación div! En cambio, hace algún tipo de magia negra con bit shifting 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 mal rendimiento. (Para x86, consulte las tablas de instrucciones de Agner Fog y la guía de microarch ).

Si conoce el divisor antes de tiempo, puede evitar la división reemplazándola por un conjunto de otras operaciones (multiplicaciones, adiciones y turnos) que tengan el efecto equivalente. Incluso si se necesitan varias operaciones, a menudo sigue siendo muchísimo más rápido que la división entera en sí misma.

Implementar el operador C de esta manera en lugar de hacerlo con una secuencia de instrucción múltiple que involucra a div es solo la forma predeterminada de GCC de dividir por constantes. No requiere optimización en todas las operaciones y no cambia nada incluso para la depuración. (Sin -Os uso de -Os para un tamaño de código pequeño hace que GCC use div .) Usar un inverso multiplicativo en lugar de división es como usar lea lugar de mul y add

Como resultado, solo tiende a ver div o idiv en la salida si no se conoce el divisor en tiempo de comstackción.

Para obtener información sobre cómo el comstackdor genera estas secuencias, así como el código que te permite generarlas para ti mismo (casi seguro innecesario a menos que estés trabajando con un comstackdor braindead), consulta libdivide .

Dividir por 5 es lo mismo que multiplicar 1/5, que es lo mismo que multiplicar por 4/5 y desplazar 2 bits por la derecha. El valor en cuestión es CCCCCCCCCCCCD en hexadecimal, que es la representación binaria de 4/5 si se coloca después de un punto hexadecimal (es decir, el binario de cuatro quintos se repite en 0.110011001100 ; consulte a continuación por qué). ¡Creo que puedes tomarlo desde aquí! Es posible que desee comprobar aritmética de punto fijo (aunque tenga en cuenta que está redondeado 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 sobre cómo funciona, explicando en términos de punto fijo. Muestra cómo funciona el algoritmo para encontrar los trabajos recíprocos y cómo manejar la división y el módulo firmado.

Consideremos por un minuto por qué 0.CCCCCCCC... (hexadecimal) o 0.110011001100... binario es 4/5. Divida la representación binaria por 4 (cambie a la derecha 2 posiciones), y obtendremos 0.001100110011... que mediante inspección trivial se puede agregar el original para obtener 0.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 , entonces 5x/4 = 1 , x=4/5 . Esto se representa entonces como CCCCCCCCCCCCD en hex para el redondeo (ya que el dígito binario más allá del último presente sería un 1 ).

En general, la multiplicación es mucho más rápida que la división. Entonces, si podemos salirse con 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 sea por una potencia de dos, pero en ese caso, por lo general, podemos simplemente convertir la división a un cambio de bit). 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 hacia cero) y luego realizamos un cambio adicional que se divide entre cuatro y trunca nuevamente Al observar el nivel de bit, está claro que podemos tratar ambos truncamientos 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 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 -2 .

(1/5) * 2 -2 <1/5 así 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 tratamos de 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 algo menos de 2 64 , terminamos con un error justo por debajo de 6/7 y después de dividir por 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, necesitamos 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, uword tiene N bits, una espada tiene 2N bits, n = numerador, 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 cambio (utilizado después de multiplicar), prec es precisión = N – e = N – shpre. El objective es optimizar el cálculo de n / d usando pre-turno, multiplicar y post-turno.

Desplácese hacia abajo hasta la figura 6.2, que define cómo se genera un multiplicador de palabras en dpa (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 se puede reducir 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 los multiplicadores de N + 1 bits en las figuras 4.1 y 4.2.

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

 ; upper 8 bytes of numerator = 2^(ℓ) = (upper part of 2^(N+ℓ)) ; lower 8 bytes of numerator for mlow = 0 ; lower 8 bytes of numerator for mhigh = 2^(N+ℓ-prec) = 2^(ℓ+shpre) = 2^(ℓ+e) numerator dq 2 dup(?) ;16 byte numerator divisor dq 1 dup(?) ; 8 byte divisor ; ... mov rcx,divisor mov rdx,0 mov rax,numerator+8 ;upper 8 bytes of numerator div rcx ;after div, rax == 1 mov rax,numerator ;lower 8 bytes of numerator 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. Observe cómo se maneja j = i / 7 (que debería ser el caso del multiplicador de N + 1).