ElBlo

Por un mundo sin divisiones

Hay una forma práctica de dividir enteros por constantes realizando sólo multiplicaciones y shifts.

El problema que queremos resolver, es dividir al número $m$ por $d$, ambos números enteros, haciendo división entera (es decir, «truncando» el resultado de la división). El número $m$ puede ser variable, pero a $d$ lo sabemos de antemano.

Este post no tiene nada novedoso. La técnica es bastante usada (por ejemplo, por compiladores), se menciona en el libro «Hacker’s Delight» y en «Division by Invariant Integers using Multiplication» está esta técnica y otras similares.

Notación: $\lfloor a \rfloor$ hace referencia a «la parte entera» del número $a$, es decir, «truncando» o descartando la parte decimal.

Recordando, el algoritmo de la división entera nos decía que para todo $m$, $d$ enteros, existen $k$ y $r$ enteros, tales que $m = d \times k + r \wedge 0 \leq r \lt d$. Es decir, podemos escribir a $m$ como un producto entre un cociente k y un divisor $d$ sumado a un resto $r$ positivo, menor a $d$.

Un ejemplo de esto sería escribir a $42$ como $4 \times 9 + 6$. Vale que $4 = \lfloor \frac{42}{9} \rfloor$.

En procesadores x86-64, usar instrucciones de división es costoso¹, y en particular, no hay ninguna que permita hacer divisiones de números enteros de forma empaquetada. Lo único que podemos hacer es dividir por potencias de 2, ya que eso es equivalente a shiftear a derecha ($m » 1 = \lfloor \frac{m}{2} \rfloor$).

Pero dejemos de pensar en computadoras y números enteros por un ratito. Supongamos que tenemos $x \in \R$, $x \gt 0$, entonces, dados $m$ y $d$ anteriores, vale que $\frac{m}{d} = \frac{m}{d}\frac{x}{x} = \frac{m\frac{x}{d}}{x}$.

Recordemos que nosotros a $d$ lo sabemos de antemano, por lo cual $\frac{x}{d}$ es algo que podemos precalcular. Si $x$ fuese divisible por $d$, entonces $\frac{x}{d}$ sería un número entero, y si $x$ además es potencia de dos, entonces dividir por $x$ es equivalente a shiftear.

Sin embargo, si $x$ es potencia de dos, va a ser divisible por $d$ si y sólo si $d$ también es potencia de dos.

Pero nos estamos acercando. Esto es relevante porque una de las divisiones en $\frac{m\frac{x}{d}}{x}$ la podríamos precalcular y la otra sería simplemente un shift a derecha.

Entonces, asumamos que $x$ es una potencia de dos, y que el resto de dividir a $x$ por $d$ es $r$, tomemos $w = d - r$. Luego vale que $x + w$ es divisible por $d$ y que $\frac{x + w}{d} = \lfloor \frac{x+w}{d} \rfloor$ es un número entero. Llamaremos $z$ a $\frac{x+w}{d}$.

Veamos cómo afecta agregar $w$ al resultado:

$$ \frac{m\frac{x+w}{d}}{x} = \frac{ m\frac{x}{d} + m\frac{w}{d}}{x} $$

Si distribuímos el $m$ entre $x + w$.

Luego:

$$ \frac{ m\frac{x}{d} + m\frac{w}{d}}{x} = \frac{m\frac{x}{d}}{x} + \frac{m\frac{w}{d}}{x} $$

El termino de la izquierda, ya habíamos concluído que era igual a $\frac{m}{d}$, y el termino de la derecha, si $x$ es lo suficientemente grande, será menor a $\frac{1}{d}$

Lo interesante viene ahora. Como sabemos que $\frac{m\frac{x}{d}}{x} = \frac{m}{d}$, el resultado entonces va a tener, a lo sumo, $\frac{d-1}{d}$ en su parte decimal². También, como vimos que $\frac{m\frac{w}{d}}{x} \lt \frac{1}{d}$, entonces vale que la parte decimal de $\frac{m\frac{x}{d}}{x} + \frac{m\frac{w}{d}}{x}$ será menor a 1.

Por lo cual, la parte entera de hacer la división con el $w$ agregado, es igual a la parte entera de hacer la división normalmente, independientemente del valor de $m$.

$$ \lfloor \frac{m\frac{x+w}{d}}{x}\rfloor = \lfloor \frac{m \times z}{x} \rfloor = \lfloor \frac{m}{d} \rfloor $$

Y eso es lo que buscamos.

Tal vez es un poco difícil de digerir todo esto, pero la idea es que dado un divisor $d$ constante, entero y positivo, podemos encontrar un $x$ potencia de $2$ y un $z$ entero positivo, tal que, para todo $m$ entero positivo y acotado, $\lfloor \frac{m \times z}{x} \rfloor = \lfloor \frac{m}{d} \rfloor$. Si sabemos que $m$ está acotado, el valor de $z$ y $x$ lo podemos precalcular. Estamos efectivamente cambiando una operación de división por una multiplicación y un shift a derecha.

Un ejemplo práctico: promediar 5 bytes.

Por ejemplo, si tenemos que tomar el promedio de 5 bytes, el valor máximo que podrían tomar es $255 \times 5 = 1275$, y querríamos dividirlo por $d = 5$. Tenemos que buscar un $x$ potencia de 2, tal que $ \frac{1275\frac{w}{5}}{x} \lt \frac{1}{5}$. Sabemos que $w$ puede tomar valores entre $0$ y $4$, por lo cual podemos acotar por $4$:

$$\frac{1275\frac{4}{5}}{x} \lt \frac{1}{5} \Rightarrow$$ $$\frac{1275\frac{4}{5}}{x}5 \lt 1 \Rightarrow$$ $$1275\times 4 \lt x \Rightarrow$$ $$5100 \lt x$$

Podemos ver que $x = 2^{13} = 8192$ cumple esta propiedad, por lo cual:

Sea $x = 2^{13}$, $w = 3$, ya que $2^{13} \equiv 2\ (\textrm{mod}\ 5)$, entonces $$ z = \frac{x + w}{d} = \frac{8192 + 3}{5} = 1639 $$

Luego

$$ \lfloor \frac{m}{5} \rfloor = \lfloor \frac{m \times 1639}{2^{13}} \rfloor = (m \times 1639)»13 $$

Dado que los números que estamos dividiendo son bastante acotados, podemos convencernos de que la división da correcta para todos los valores con un simple programa en Python:

>>> all([x//5 == ((x*1639)>>13) for x in range((255*5)+1)])
True

En resumen, para dividir un número entero $m$ en el rango $[0, \textrm{max}]$ por una constante $d$, entera y positiva, tenemos que:

1. Encontrar un $x$ potencia de $2$, tal que $\frac{\textrm{max} \times (d-1)}{x} \lt \frac{1}{d}$
2. Definir $n$ como lo que le falta a $x$ para ser divisible por $d$.
3. Calcular $z = \frac{x+n}{d}$

Luego, la división es simplemente hacer $\lfloor \frac{m}{d} \rfloor = \lfloor \frac{m \times z}{x} \rfloor$

Cómo usar esto en el contexto de SIMD

Si nuestros números están en words, podemos usar la instrucción pmulhw para multiplicar nuestro número por $z$. Quedándonos con la parte alta de la multiplicación, lo cual sería equivalente a hacer $(x\times z)»16$

Para esto, necesitariamos entonces elegir como $m = 2^{16} = 65536$, entonces $n = 4$, $z = \frac{65540}{5} = 13108$. El código en assembler quedaría algo así:

section .rodata
ALIGN 16
div5: times 8 dw 13108

section .text
; dividir_por_5 divide 8 words sin signo en xmm0 por 5.
; Los words tienen que estar entre 0 y 255*5.
dividir_por_5:
  pmulhw xmm0, [div5]
  ret