The Quake3’s Fast InvSqrt()

1.Introdução

O inverso da raiz quadrada é uma função presente em diversos problemas matemáticos, físicos e estatísticos. Aplicações algébricas de normalização, desvio padrão e até na indústria de videogames possuem freqüente utilização da simples expressão:

E foi na programação de um jogo que surgiu uma função que acabou por se tornar famosa entre os programadores por conseguir rodar até 4 vezes mais rápido que a implementação derivada da função sqrt(x), nativa da biblioteca do C.

A função que realizava o cálculo tão mais rapidamente com um erro associado de 1,75x10E-4 se apresentou muito eficiente na época e significou uma melhora considerável na execução do Quake – jogo que utilizava esta operação frenquentemente e para qual foi desenvolvida. O seu código é:

float InvSqrt(float x){
	float xhalf = 0.5f*x;
	int i = *(int*)&x;

	i = 0x5F3759DF – (i>>1);
	x = *(float*)&i;
	x = x*(1.5f – xhalf*x*x);

	return x;
}

Sua verdadeira autoria é desconhecida. Contudo, tudo indica que foi programada por John Carmack, na época programador do Quake 3, onde surgiu.

Para entender esta a genial “sacada” é necessário conhecer o algoritmo matemático usado e como ele foi optimizado graças ao conhecimento do padrão de ponto flutuante pelo autor.

2.O Método de Newton-Rapson

O algoritmo utilizado foi um método de aproximações numéricas já bastante conhecido em muitos projetos computacionais, conhecido como “Aproximação de Newton”, ou “Método de Newton-Rapson”.

Este é um método linear iterativo utilizado para encontrar raízes x de uma função f(x) a partir de uma estimativa inicial baseado na idéia de usar retas tangentes para substituir o gráfico y = f(x) próximo de onde f é zero.

Para tal, cria-se uma função G, chamada “função iteração” xnew = G( xold ) , a partir da função f(x) para se chegar à raiz de f(x). Sendo assim, a convergência deste método será tão rápida quanto for |G'(xold)|.

O método de Newton-Rapson modifica o método iterativo linear simples (Método do Ponto Fixo Modificado) fazendo com que |G'( xold )|->0. Assim, a fórmula iterativa fica:

Uma propriedade importante para a função InvSqrt é que o Método de Newton-Rapson converge quadraticamente. Uma interpretação gráfica do método pode ser vista na figura abaixo:

2.1.Newton-Rapson na InvSqrt

Queremos o retorno da função y=y(x) tal que:

contudo, queremos o valor de y, que é a inversa da raíz quadrada. Sendo assim, procuramos a função inversa tal que:

como sabemos que o Método de Newton-Rapson para encontrar raízes f(x) = 0, então poderíamos substituir a variável x pela y, que é a que queremos encontrar.

como

então

aplicando o Método de Newton, teremos que:

onde sabemos f(y),

logo

finalmente, a expressão para encontrar o valor de y é:

que pode ser reduzida para:

ou passando 0.5 para expressão dentro dos parentes,

lembrando que quanto maior o número de iterações, mais próximo da raiz será o valor do resultado.

3.I3E 754 de 1985

Para entender a maior “mágica” deste código, é necessário saber como os números de ponto flutuante são armazenados na memória e interpretados pelo processador.
Em ciências exatas é comum se ver a notação científica de base 10 para expressar valores muito grandes ou muito pequenos. Esta notação consiste em usar apenas um único dígito à esquerda do valor decimal. Por exemplo, um valor em notação científica na base 10 é:

0,8 x 10^-8

Da mesma forma, podemos expressar números em notação científica em base 2. Por exemplo:

1,0 x 2^7 (em base 2)

Contudo, como não se está trabalhando com a base 10, é conveniente renomear os valores após a vírgula de decimais para ponto binário.
A expressão “ponto flutuante” vêm porque a representação dos números em que o ponto binário não é fixo – por isso flutua. A linguagem em que o InvSqrt foi implementada, C, utiliza esta o nome “float” para estes números de ponto flutuantes.
Desta forma, internamente o computador armazena todos expoentes em base 2, algo como:

1,0 x 10^111 ( binário equivalente à 1,0 x 2^7 )

Observe que a representação de ponto flutuante implica em encontrar uma relação entre o tamanho da fração e o tamanho do expoente, pois o número de bits da palavra do processador é fixo. Assim, quando se aumenta o tamanho da fração, reduz-se o valor do expoente.
Logo, os números em ponto flutuante seguem o formato:

(-1)^s x F x 2^E

onde s é o sinal do número de ponto flutuante de maior relevância e que também é o bit que indicará o sinal do número (s=0 implica em número positivo, s=1, negativo).

F envolve o campo de fração e E são os bits de expoente. Com base nessa lógica, foi elaborado um padrão para pontos flutuantes chamado “padrão IEEE 754”, que é implementado em quase todos processadores fabricados desde 1985. A adoção deste padrão permitiu uma portabilidade muito maior entre plataformas e melhorou muito a qualidade numérica de operações aritméticas com valores reais.
Para entender o padrão, veja na tabela abaixo como se relacionam os bits na palavra de uma arquitetura de 32 bits:

Note que reservando 8 bits para o expoente temos que 7 bits são utilizados para o valor do expoente e 1 sinal é reservado para o sinal do expoente. Isso equivale em decimal à um range de 2×10^-38 a 2×10^38.
Seguindo a mesma idéia da notação de ponto flutuante de precisão simples (float), é implementada a precisão dupla (tipo “double”). Esta porém, utiliza 2 registradores de 32 bits em arquiteturas com esta palavra de processador, ou 1 registrador de 64 bits para a representação de ponto flutuante com precisão maior.
No caso da precisão dupla em processadores 32 bits, 1 registrador armazena 32 bits para a fração e outro registrador possui a seguinte divisão:

Ou seja, tanto em arquiteturas de 64 quanto de 32 bits são reservados 52 bits para fração (20 mesmo registrador + 32 de um segundo registrador em processadores 32 bits), 11 bits para expoente (1 de sinal + 10 de valor) e 1 bit para o sinal do número.
Logo, a precisão dupla permite um range de valores de 2×10^-308 a 2×10^308.
Outro ponto importante para compreensão do InvSqrt é saber que os projetistas do IEEE 754 construíram esta representação para que se pudesse ser facilmente processada utilizando os mesmos circuitos de comparação de inteiros – provavelmente já implementados nas arquiteturas da época. O motivo disso está relacionado com o melhor aproveitamento dos registradores (que não precisam ser dedicados à valores de ponto flutuante ou valores inteiros) e para permitir um teste rápido de comparação (instruções bne e beq aproveitadas para qualquer tipo de dado).

Assim, a representação do IEEE 754 pode ser processada com as mesmas comparações de inteiros para acelerar a ordenação dos números de ponto flutuante.
O IEEE 754 é sem dúvida um padrão muito consistente e eficiente e vêm sendo mantido por sua confiabilidade. Além das propriedades descritas aqui, ainda há algumas outras que vale a pena conhecer. Por exemplo, o padrão já prevê um símbolo padrão para operações que geram erros, utilizado por “baixo dos panos” no envio de sinal para o programa tratar uma exceção. Uma leitura altamente recomendada deste padrão está no livro “Computer Organization and Design”, do Patterson.

4.InvSqrt

Começando pela primeira linha do algoritmo, pode-se adiantar que xhalf equivale ao produto (0.5*x) da expressão obtida no ítem 2.1:

float xhalf = 0.5f*x;// equivalente a expressao acima

este valor é salvo em uma nova variável (xhalf) porque a variável x será reaproveitada para guardar o valor de raiz (chamada de y nas demonstrações). Como visto anteriormente, a variável x (das expressões) só aparece na última expressão (0.5*x).
Na linha seguinte:

int i = *(int*)&x;

temos um casting de ponteiros que é a alma deste código. Internamente, o que ocorre é uma mudança na significância dos bits.
Enquanto x é uma variável de ponto flutuante de 32 bits, com a significância para os bits do padrão 754 (1 bit de sinal, 8 para o expoente e 23 para a base), a variável i é uma variável inteira, onde todos os bits possuem a mesma significância. Todavia, como tanto o tipo de x quanto de i possuem o mesmo tamanho em palavra – 32 bits – é possível armazenar ambas em um mesmo espaço. Isso é o que faz esse código ser possível.
As operações com ponteiros utilizada aqui deve ser lida de trás pra frente. A ordem diz que:

– &x: passe o endereço do espaço de memória onde x ocupa (4 palavras de memoria, ou 32 bits), que suporta tipos de ponto flutuante.

– (int *): converta este espaço de memória para suportar tipos de dados inteiros simples (também 32 bits)

– *: retorne o valor armazenado lá como inteiro.

Na próxima linha,

i = 0x5f3759df – (i>>1);

têm-se a estimativa inicial para y[n] = y[0]. Note que a operação (i >> 1), deslocamento de bits à direita, equivale à divisão por 2. O uso desta operação é a justificativa para a conversão da representação de float para int.
Veja que realizando o deslocamento de bits em todo o vetor de bits, também se tira a raíz do expoente. Por exemplo, o deslocamento da representação por inteiro de um valor qualquer:

BIN(01100110101111100010110101010101) = DEC(3447478954)

apenas deslocando os bits à direita 1 vez, ( 01100110101111100010110101010101 >> 1) têm-se que:

BIN(00110011010111110001011010101010) = DEC(1723739477)

onde se vê claramente que 1723739477 = 3447478954/2.
Mas isso não tem nada de excepcional. Contudo, analisando levando em conta a reserva de bits para sinal, base e expoente, da norma IEEE 754 para essa mesma seqüência, vê-se no exemplo:

quando é aplicada o deslocamento de bits à direita, têm-se:

observe que o expoente é exatamente a metade inteira do valor antes do deslocamento dos bits. Da segunda linha das tabelas, no valor decimal, $$(int) 102.5 == 102 == (int)205/2 $$. Portanto, o motivo do autor do código de converter os bits de uma variável float para uma inteira e deslocar o bits é justamente para dividir o expoente por 2, que obviamente equivale à tirar a raiz. Em seguida, quando o autor reconverte o valor para ponto flutuante,

x = *(float*)&i;

ele simplesmente recoloca os bits na forma interpretada pelo processador como sendo um float do IEEE 754.

Sendo assim, após já tirar a raíz do valor de entrada, o método de Newton é aplicado na última linha antes do retorno da função, para convergir ao valor da raiz quadrada:

x = x * (1.5f – xhalf * x * x);
// xhalf = 0.5x e x = y, nas demonstracoes

conforme demonstrado, o Método de Newton para raiz quadrada inversa:

5.Conclusões

O código da InvSqrt quando foi desenvolvido, rodava supostamente 4 vezes mais rápido que utilizando (float) 1/sqrt(x), da biblioteca padrão e o erro relativo está por volta de 10^-4. O que tornava o código muito útil.
Com o avanço das otimizações dos compiladores e dos circuitos de multiplicação dos processadores, a função nativa consegue ser mais eficiente. Contudo, conhecer esta função, seu funcionamento e o algoritmo utilizado é engrandecedor para qualquer programador.
A importância desta função foi tamanha que muita discussão foi gerada em torno dela. Inclusive modificações, melhorias e até alguns artigos sobre ela.
Um exemplo disso foi a pesquisa de Chris Lomont (www.math.purdue.edu/~clomont) sobre o melhor valor a ser usado como y[0] no número mágico para que a primeira aproximação tenha convergência mais rápida. Segundo ele, a melhor constante para melhoria da performance seria 0x5f375a86, uma pequena diferença de 167 em relação ao original.

6.Referencias

• Finney, Ross L. CÁLCULO DE GEORGE B. THOMAS JR. , volume 1. Tradução de Paulo Buschcov. Addison Wesley, 2002. Páginas 301-305. “O MÉTODO DE NEWTON”.
• Chris Lomont, www.math.purdue.edu/~clomont, “FAST INVERT SQUARE ROOT”.
• IEEE 754, grouper.ieee.org/groups/754
• Patterson, David A. COMPUTER ORGANIZATION AND DESIGN. 3rd ed. Elsevier, 2005.