Revisitando la raíz cuadrada invertida

Hace algunos años se puso de moda hablar en blogs del algoritmo que se incluyó en el Quake III Arena para calcular una aproximación rápida de la raíz cuadrada invertida. Esta operación matemática es útil para calcular los ángulos de incidencia y reflexión de los efectos de iluminación y profundidad.

De hecho se ha hablado ya mucho sobre el tema así que para ahorraros tiempo quiero comentar que hay una instrucción máquina en el conjunto SSE (Streaming SIMD Extensions) que la calcula mucho más deprisa que si escribimos nosotros el código. Dicha instrucción es la rsqrt de hecho si no buscamos una raíz cuadrada exacta calcular rsqrt(x) * x puede llegar a ser más rápido que la instrucción sqrt, aunque esta última tiene más precisión. De todas formas no está de más probar las cosas siempre para estar totalmente seguros.

Si queréis profundizar en el tema. el algoritmo se discute en este artículo de investigación escrito en el año 2003.

float fastInvSqrt(float x) {
    float xhalf = 0.5 * x;
    int i = *(int *)&x; // Castea a entero el número en coma flotante
    i = 0x5f3759df - (i >> 1); // Número mágico
    x = *(float*)&i; // Vuelve a castear a float
    x = x*(1.5-(xhalf*x*x)); // Aplica una iteración del método de Newton
    return x;
}

El número mágico permite aproximar la raíz cuadrada inversa usando aritmética entera, para obtenerlo se usan propiedades de logaritmos y cambios de representación entre enteros y coma flotante. Después se aplica una iteración del método de Newton.

Podéis encontrar una descripción completa en este enlace de wikipedia.

Artículos relacionados:

Contando bits

Guy Kawasaki, evangelismo digital

Allá por el año 1992 llegó mi primer ordenador a casa, un Macintosh Performa 460. ¿Como llegó este ordenador a mi habitación? Bueno en aquella época Linux acababa de nacer, no estaba extendido. Además, las versiones de Windows anteriores a Windows 95 eran muy inferiores a lo que el MacOS 7.5 podía ofrecer. Por último, el software / hardware gráfico y de sonido era el mejor (Apple Sound Manager System), disfrutando de las mejores versiones de juegos como Sim City 2000 o Prince of Persia (flame on).

Pero ninguna de estas razones fueron el motivo por el que este ordenador entró en casa, ya que yo tenía 12 años y no sabía nada de informática y mis padres menos. La culpa de todo la tiene este señor Guy Kawasaki, nacido en Hawai en 1954:

Guy Kawasaki. Licencia: CC by 2.0. Fuente: Wikipedia, GeeJo

Guy Kawasaki entró a trabajar en 1983 en Apple Computer. Allí trabajó como evangelista en jefe durante 4 años, haciendo marketing del Macintosh. Volvería después a trabajar para Apple en 1995. El término evangelista lo acuñó Mike Murray, otro trabajador de Apple. Mediante el evangelismo convenció a mi vecino de que ese era el ordenador que debíamos comprar.

¿En qué consistía el evangelismo? Pues bien era una técnica de marketing en la que se convence y entusiasma a los usuarios de tal manera con un producto que ellos mismos lo van a recomendar a otras personas. Hoy en día, esto se conoce también como canales de confianza (trust), dentro del marketing talk, think and trust. Pero aquello tenía un rollo más espiritual, como una religión. Las revistas del sector y las páginas de internet (aún no existían los blogs) seguían este modelo en sus contenidos, fidelizando a los lectores hacia la marca y posicionándolos en contra del PC. Además lo reforzaban con anuncios como el de Think different y ya te creías de verdad que estabas salvando a tus colegas al recomendarles un Mac:

La verdad fue bonito mientras duró, y muy eficaz. Reconozco ser una víctima de aquella estrategia y que te mandaría a comprar un Mac de por entonces sin pensarlo. Pero salió Windows 95, y llegó el software libre y aquello terminó. Personalmente pienso que fue un arma de doble filo y al crear un vínculo tan personal muchos evangelistas se sintieron abandonados por Apple. Incluso ediciones españolas de revistas como MacFormat, muy centradas en la evangelización, tuvieron que cerrar. Recordemos que a la firma le fue muy mal hasta que Steve Jobs volvió, ya que sus ordenadores perdieron la ventaja de la innovación.

Sobre esto Guy Kawasaki dijo en una entrevista que si alguna vez llegaba a Apple otra persona como Steve pero con las ideas incorrectas sobre los productos a fabricar hundiría a la compañía. Además comenta en sus charlas que es un milagro que Apple siga con vida.

Pero lo que más nos importa es que Guy Kawasaki supo evangelizar y vender, y que ahora da charlas muy buenas y escribe libros para emprendedores:

Lo que más destaca es su estilo transgresor para introducir conceptos de marketing. Por ejemplo: "Don't worry be crappy" (No te preocupes se mierdoso), con la que exhorta a sacar el producto rápido para aprovenchar la ventaja competitiva.

¿Qué os parece? ¿Habéis sido alguna vez víctimas del evangelismo? ¿Habés tenido a algún amigo pesado, víctima del marketing?.

Contando bits

Las operaciones de cuenta de bits tienen muchas aplicaciones, como por ejemplo compresión de datos. En el ordenador todos los valores numéricos se codifican en sistema binario. Así pues, usando 32 bits de precisión el número 212 se codificaría como:

0000 0000 0000 0000 0000 0000 1101 0100

en un registro del procesador. La operación que cuenta el número de bits a 1 en un registro se conoce como bitcount() o popcount(). También se puede aplicar la definición del peso de Hamming, que es el número de elementos de una cadena distintos de 0. En el caso del ejemplo:

bitcount(212) = 4

La implementación de esta operación depende de la arquitectura de la CPU que estemos utilizando. En ocasiones, la CPU tendrá una instrucción hardware dedicada al bitcount(). Las instrucciones POPCNT y LZCNT forman parte del set ABM (Advanced Bit Manipulation). Están incluidas en la arquitectura x86 en el conjunto de instrucciones SSE 4.2.

Tener una instrucción hardware dedicada es el caso más óptimo, ya que en pocos ciclos de procesamiento tendremos nuestro resultado. Si estamos programando en un lenguaje de alto nivel debemos buscar los pasos para usar la instrucción máquina en nuestra arquitectura (normalmente llamamos a una función específica y compilamos con un flag especial).

Sin embargo existen casos en los que no disponemos de esta instrucción. En estos casos una buena solución general consiste en usar operaciones bit a bit (esta entrada de la wikipedia es fenomenal para entenderlo y antes de seguir también deberíais dominar la representación de bits en formato hexadecimal). Si estamos trabajando con registros sin signo de 32 bits el código sería el siguiente:

inline uint32_t popcount(uint32_t i) {
    i = i - ((i >> 1) & 0x55555555);
    i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
    return (((i + (i >> 4)) & 0xF0F0F0F) * 0x1010101) >> 24;
}

Este algoritmo se conoce como HAKMEM. La culpa de todo es de unos tales Beeler, Gosper and Schroppel, allá por el año 1972. Como véis la función es un poco confusa, así que podéis añadir un comentario diciendo "// Funciona" o "// Ojo con esta pomada". O podéis seguir leyendo la entrada y crujirme en los comentarios si no lo explico bien.

HAKMEM emplea una estrategia divide y vencerás para contar en paralelo. Veamos como lo hace tomando como ejemplo el número 212.

1) Desplazamos un bit a la derecha el registro inicial i >> 1:

0000 0000 0000 0000 0000 0000 1101 0100
----------------------------------------
0000 0000 0000 0000 0000 0000 0110 1010

2) Ahora aplicamos la función AND con el valor 0x55555555 ( i >> 1 ) & 0x55555555:

0000 0000 0000 0000 0000 0000 0110 1010

0101 0101 0101 0101 0101 0101 0101 0101
----------------------------------------
0000 0000 0000 0000 0000 0000 0100 0000

3) El resultado se lo restamos al valor inicial i - ( ( i >> 1 ) & 0x55555555 ):

0000 0000 0000 0000 0000 0000 1101 0100

0000 0000 0000 0000 0000 0000 0100 0000
----------------------------------------
0000 0000 0000 0000 0000 0000 1001 0100

Con estas tres operaciones iniciales lo que hemos obtenido es el número total de 1 en cada par de bits del registro original. Si nos fijamos en los últimos 8 bits del registro original y nuestro resultado actual:

11 01  01 00

10 01  01 00

Podemos ver que el par 11 tiene 2 bits a uno (10) y que le siguen dos pares 01 con 1 bit a uno (01) y un par 00 con ningún bit a uno (00).

La siguiente operación se encarga de sumar estos valores de dos en dos, de forma que vamos a obtener el número de bits a uno en cada grupo de cuatro bits, veamos como:

4) Primero aplicamos la máscara 0x33333333 a nuestro último valor i & 0x33333333:

0000 0000 0000 0000 0000 0000 1001 0100

0011 0011 0011 0011 0011 0011 0011 0011

----------------------------------------

0000 0000 0000 0000 0000 0000 0001 0000

5) Desplazamos el último valor dos bits a la derecha y aplicamos la máscara (i >> 2) & 0x33333333:

0000 0000 0000 0000 0000 0000 1001 0100
----------------------------------------
0000 0000 0000 0000 0000 0000 0010 0101

0011 0011 0011 0011 0011 0011 0011 0011

----------------------------------------

0000 0000 0000 0000 0000 0000 0010 0001

6) Sumamos estos dos valores:

0000 0000 0000 0000 0000 0000 0001 0000

0000 0000 0000 0000 0000 0000 0010 0001

----------------------------------------

0000 0000 0000 0000 0000 0000 0011 0001

Ya tenemos el número total de unos en cada grupo de cuatro bits. Si comparamos con el valor original de los últimos 8 bits del registro:

1101 0100

0011 0001

Podemos ver que el primer grupo 1101 tiene 3 bits a uno (0011) y el segundo grupo 0100 tiene un bit a uno (0001).

El siguiente grupo de operaciones obtiene el número total de unos en cada grupo de 8 bits.

7) Sumamos i + (i >> 4):

0000 0000 0000 0000 0000 0000 0011 0001

0000 0000 0000 0000 0000 0000 0000 0011

----------------------------------------

0000 0000 0000 0000 0000 0000 0011 0100

8) Aplicamos la máscara de bits 0xF0F0F0F (i + (i >> 4)) & 0xF0F0F0F:

0000 0000 0000 0000 0000 0000 0011 0100

0000 1111 0000 1111 0000 1111 0000 1111

----------------------------------------

0000 0000 0000 0000 0000 0000 0000 0100

Si comparamos con el valor original de los últimos 8 bits (byte) del registro:

11010100
00000100

Vemos que concuerda que 11010100 contiene 4 bits (00000100).

Creo que ya os podéis hacer una idea de como funciona el algoritmo. Para nuestro ejemplo ya habríamos terminado, porque sólo usa los 8 bits menos significativos. Sin embargo, si habéis tenido la paciencia de llegar hasta aquí os vais a llevar una gran sorpresa en forma de propiedad matemática que nos permite obtener la cuenta de todo el registro de 32 bits.

Ahora mismo cada grupo de 8 bits (un byte) contiene la cuenta de bits a uno en el registro original para ese byte. Multiplicar por 0x1010101 tiene una propiedad especial, y es que si el número por el que multiplicamos tiene 4 bytes (A, B, C, D) el resultado de la multiplicación nos devolverá un resultado cuyos 4 bytes serán (A+B+C+D, B+C+D, C+D, D). Es decir, el primer byte contiene la cuenta de bits del registro de 32 bits, por lo que ahora sólo nos quedaría desplazar el contenido de este byte 24 bits a la derecha >> 24 para obtener el resultado.

Y bien, eso ha sido todo por hoy, espero que os haya gustado esta entrada tanto como me ha gustado a mi prepararla. Comentarios, sugerencias, ríos de tinta... son bienvenidos.

Algunas referencias:

Stack overflow

Página de la universidad de Standford

Hackers delight, capítulo 5

Operaciones bit a bit

Entradas relacionadas:

Revisitando la raíz cuadrada invertida

Hola a todos

Bienvenidos a mi blog, donde junto al androide Chipperman vamos a disfrutar al máximo de la informática y las ciencias de la computación. Chipper en inglés significa alegre, contento y animado, así es como me siento yo empezando este proyecto del que vosotros con vuestros comentarios y aportaciones vais a ser la pieza clave. Tocaremos temas de gran interés, tanto en el apartado técnico como en el de la organización y consecución de proyectos.

Lo que hace a Chipperman un androide perfecto no es sólo su actitud y gran motivación, Chipper en inglés también significa astilladora. Él es capaz de coger esos problemas que a veces no entendemos y triturarlos. Será nuestra tarea coger ese serrín y darle forma. Y es que lo que hace grande a un informático no es dominar muchos lenguajes de programación con el máximo detalle, sino entender los distintos problemas abordados en computación y las herramientas matemáticas y algoritmos que se utilizan para resolver cada uno de ellos.

Perfect Chipperman es también el nombre de una canción de un músico de la escena de Atari ST, MC Laser / Lotek Style. Aquí tenéis los enlaces a la canción y a la página del autor. Necesitaréis un reproductor de ficheros SNDH para escucharla.

La cosa irá adelante, de momento ya tengo un canal Youtube donde voy a ir subiendo tutoriales.