- O fast inverse square root, que ficou famoso com o Quake 3, era uma solução de desempenho da época para aproximar rapidamente
1 / sqrt(x) por meio da reinterpretação dos bits de um float e de uma correção de Newton-Raphson
- O ponto central é que o padrão de bits inteiros de um float IEEE-754 de 32 bits pode ser tratado como uma aproximação escalada e deslocada de
log2(x)
0x5f3759df - (i >> 1) é uma forma de traduzir log2(x^-0.5) = -0.5 * log2(x) para deslocamento e subtração inteiros, e a constante mágica vem de 3/2 * 2^23 * (127 - σ)- Em seguida,
y = y * (1.5 - 0.5x * y * y) aplica uma correção de Newton-Raphson uma vez, e a segunda iteração no código do Quake está comentada
- Em 1999, a inversa da raiz quadrada era necessária centenas ou milhares de vezes por segundo para iluminação e normalização de vetores 3D, mas em hardware moderno a utilidade prática desse mesmo truque diminuiu graças ao processamento dedicado de ponto flutuante
O que o código do Quake faz
float Q_rsqrt(float number) {
long i;
float x2, y;
const float threehalfs = 1.5F;
x2 = number * 0.5F;
y = number;
i = *(long*)&y;
i = 0x5f3759df - ( i >> 1 );
y = *(float*)&i;
y = y * ( threehalfs - ( x2 * y * y ) );
return y;
}
- Essa função calcula um valor aproximado da inversa da raiz quadrada
1 / sqrt(number) para number
- A parte mais famosa é a manipulação de bits que interpreta o valor float como
long e depois executa 0x5f3759df - (i >> 1)
- Quando o Quake 3 foi lançado, em 1999, a inversa da raiz quadrada era uma operação lenta e cara, mas era necessária centenas ou milhares de vezes por segundo em cálculos de vetores 3D que exigiam equações de iluminação e normalização
- Em hardware moderno, esse tipo de cálculo pode nem ser executado na CPU ou, mesmo quando é, acaba sendo rápido graças ao hardware dedicado de ponto flutuante, muito mais evoluído
Representação IEEE-754 de float de 32 bits
- Um float de 32 bits é composto por três partes
- Sign: 1 bit, indica se o valor é positivo ou negativo
- Exponent: 8 bits, define a faixa em que o valor se encontra
- Mantissa: 23 bits, representa linearmente a posição dentro dessa faixa
- Um valor comum é interpretado na forma abaixo
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
B = 127 é o valor de viés usado no biased exponent, e o expoente real é e = E - B- A mantissa não é usada simplesmente como multiplicação por
m, mas na forma 1 + m
- se
m = 0, temos 2^e
- quando
m se aproxima de 1, representa valores até pouco antes da próxima faixa de expoente, 2^(e+1)
- Se todos os bits do expoente forem 0, trata-se de um número sub-normal, e a fórmula muda
N = (-1)^S * 2^-126 * m
- Os sub-normais são necessários para representar 0 e números muito pequenos, próximos de 0
- Se todos os bits do expoente forem 1, o valor é tratado como especial
- se
E = 255 e M = 0, temos Infinity ou -Infinity
- se
M != 0, temos NaN
A relação logarítmica que surge ao ver os bits do float como inteiro
- Se a representação interna do float for vista como um inteiro de 32 bits, ela pode ser escrita assim
I_x = 2^31 S + 2^23 E + M
- Como a inversa da raiz quadrada é aplicada a entradas positivas, podemos considerar
S = 0, o que simplifica a expressão
L = 2^23
I_x = L E + M
- Dentro da mesma faixa de expoente, a mantissa representa a posição de forma linear, mas conforme o expoente cresce, a mesma quantidade de passos da mantissa cobre um trecho maior da reta numérica
E = 127, ou seja, e = 0, corresponde aproximadamente ao intervalo [1, 2)E = 128, ou seja, e = 1, corresponde aproximadamente ao intervalo [2, 4)- os dois intervalos têm o mesmo número de passos de mantissa, mas o segundo é duas vezes mais largo
- Por causa dessa estrutura, quando o padrão bruto de bits do float é visto como inteiro, surge uma relação logarítmica
Os bits brutos como aproximação de log2(x)
- Se o padrão de bits de um float for interpretado como o inteiro
I_x, ele pode ser visto como uma aproximação linear por trechos de log2(x)
- Essa relação pode ser expressa pela aproximação abaixo
log2(x) ≈ I_x / L - B
- Se dividirmos o inteiro bruto dos bits pelo tamanho da mantissa
L = 2^23 e subtrairmos o viés do expoente B = 127, obtemos um valor próximo de log2(x)
- O logaritmo dentro do intervalo da mantissa é tratado por aproximação linear
log2(1 + x) ≈ x + σ
σ é um parâmetro de ajuste da aproximação, e x representa a posição dentro da faixa de expoente, no intervalo [0, 1]
Transformando a inversa da raiz quadrada com uma identidade logarítmica
- O objetivo é calcular o valor abaixo
y = 1 / sqrt(x)
- Isso pode ser reescrito na forma exponencial a seguir
y = x^-0.5
- Aplicando a identidade logarítmica, o cálculo da inversa da raiz quadrada vira a relação abaixo
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
- Usando o fato de que os bits do float se comportam como uma aproximação de
log2(x), podemos aproximar diretamente a representação inteira I_y de y a partir da representação inteira I_x de x
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
- Essa expressão leva diretamente à linha central do código do Quake
i = 0x5f3759df - ( i >> 1 );
i >> 1 desloca os bits inteiros uma posição à direita e cumpre o papel de multiplicar por 1/2- A constante à frente,
0x5f3759df, corresponde a 1.5 * L * (B - σ)
O que é a constante 0x5f3759df
- Se fizermos
σ = 0, a constante é calculada da seguinte forma
1.5 * 2^23 * 127 = 1598029824
- A representação hexadecimal desse valor é
0x5f400000
- Há uma diferença de
566817 em relação à constante real usada no Quake, 0x5f3759df
- A partir dessa diferença, podemos calcular o valor de
σ correspondente ao código do Quake
σ = 377878 / 2^23
σ = 0.04504656
- Em C, a mesma constante pode ser calculada assim
int32_t compute_magic(void) {
double sigma = 0.0450465;
double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
int32_t i = expression;
return i;
}
// -> 0x5f3759df
- Aqui, usa-se
double, e a conversão para inteiro não é uma reinterpretação de bits, mas um cast comum
- Esse valor de
σ foi escolhido para otimizar a aproximação, mas não é o valor realmente ótimo, e também não se sabe com certeza quem o criou
Por que isso não é só um hack simples
0x5f3759df - (i >> 1) é uma expressão que cria um valor inicial para a inversa da raiz quadrada explorando o fato de que os bits brutos do float funcionam como aproximação de logaritmo- Ela se apoia em uma relação matemática complexa, mas na execução usa apenas operações rápidas como deslocamento e subtração
- Na época, era necessário processar operações caras milhares de vezes por segundo, então essa abordagem virou um projeto de engenharia adequado às limitações do hardware
- No entanto, esse algoritmo só funciona com normal float
- em valores sub-normais, a hipótese da aproximação
log2(1 + x) ≈ x + σ não vale
- em sub-normais, o formato real fica mais próximo de
0 + x, o que quebra a aproximação
Reduzindo o erro com correção de Newton-Raphson
- O valor inicial obtido pela manipulação de bits já é muito bom, mas ainda deixa um erro mensurável
- A linha abaixo melhora bastante a aproximação
y = y * ( threehalfs - ( x2 * y * y ) );
- Essa linha aplica o método de Newton-Raphson
- Para adaptar o problema da inversa da raiz quadrada ao método de Newton, ele é transformado em um problema de encontrar a raiz da função abaixo
f(y) = 1 / y^2 - x = 0
- O método de Newton gera, a partir da aproximação atual
y_n, uma aproximação melhor y_(n+1) da seguinte forma
y_(n+1) = y_n - f(y_n) / f'(y_n)
- A derivada de
f(y) = y^-2 - x é a seguinte
f'(y) = -2y^-3 = -2 / y^3
Fórmula de correção de Newton sem divisão
- Se aplicarmos a fórmula de Newton diretamente, várias divisões em ponto flutuante aparecerão
- Um dos motivos de esse algoritmo ser rápido é justamente a evasão de divisões em ponto flutuante
- Fazendo a simplificação algébrica, chegamos a uma forma que usa apenas multiplicações, sem divisão
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
- No código do Quake,
x2 = number * 0.5F calcula 0.5x de antemão, e isso é usado na linha seguinte
y = y * ( threehalfs - ( x2 * y * y ) );
- Após essa única iteração, o erro absoluto máximo é de 0,175%, e em muitos casos o erro fica abaixo disso
- O código original inclui uma segunda iteração de Newton, mas ela está comentada
// y = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed
Origem e algoritmos relacionados
- Esse algoritmo não foi inventado por John Carmack, e sua origem exata não é 100% certa
- Há um texto relacionado da Beyond3D com este link: The truth is the exact origin is not 100% certain
- Chris Lomont escreveu um artigo tentando encontrar o valor ótimo de sigma na etapa de aproximação logarítmica: InvSqrt.pdf
- CORDIC é um algoritmo que calcula seno e cosseno usando apenas somas e deslocamentos de bits, sem ponto flutuante, e difere bastante do fast inverse square root nos detalhes
- Os dois algoritmos têm em comum o fato de aplicarem observações matemáticas de maneira eficiente, de acordo com as limitações de hardware da época
2 comentários
Aquele código curioso que reaparece de tempos em tempos... haha
Comentários no Hacker News
Se o computador foi fabricado depois de 1999, em geral ele oferece suporte ao conjunto de instruções SSE, que inclui
_mm_rsqrt_ps, capaz de calcular quatro raízes quadradas inversas de uma vez mais rapidamente: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...Ainda assim, a técnica discutida aqui ainda não é completamente inútil. Há hardware em que conversões float/int são rápidas, mas que não tem instruções
rsqrt,sqrt,powoulog, e essas operações podem ser aproximadas com esse truquesqrt(x)comox * 1/sqrt(x)Em conjuntos de instruções de GPU, ARM, RISC-V, AVR, PIC, 8051, FPGA etc., muitas vezes há uma operação aproximada de raiz quadrada inversa embutida, mas é bem provável que ela seja implementada com algoritmos desse tipo
Para implicar um pouco com o artigo: a explicação de que esse tipo de cálculo não acontece em CPUs modernas não está correta. É um equívoco comum achar que jogos ou apps com muitas operações de ponto flutuante querem mandar todas as operações de ponto flutuante para a GPU
Na prática, só faz sentido mandar para a GPU trabalhos grandes e uniformes. Se você está fazendo uma normalização pontual de vetor, como montar uma matriz de rotação para fazer um objeto olhar para outro, é mais rápido deixar isso na CPU. Mesmo desconsiderando o tempo de transferência para a GPU, uma única operação de ponto flutuante é mais rápida na CPU, porque GPUs normalmente têm clocks mais baixos e alcançam FLOPS altos por meio de paralelismo
Escrevi uma implementação em MMIX e parti da suposição de que o valor de entrada original é maior que
2^-1021Se houver interesse, a Wikipedia também tem uma boa explicação sobre essa função e sua história: https://en.wikipedia.org/wiki/Fast_inverse_square_root
Reuni algumas coisas desse tipo aqui: https://github.com/ncruces/fastmath/blob/main/fast.go
Há também uma discussão relacionada no StackOverflow: https://stackoverflow.com/questions/32042673/optimized-low-a...
fastmathChegou a hora de implicar. Há um erro de digitação na fórmula de float: deve ser
(-1)^S, não-1^S. Esta última é sempre-1A explicação de que interpretar o padrão bruto de bits é uma aproximação linear por partes do logaritmo também não é precisa. As linhas entre os pontos de dados no gráfico azul não existem de fato, e não dá para ter um bit definido como 1 pela metade. É mais parecido com uma versão discreta do logaritmo, e os pontos de dados que realmente existem — isto é, os pontos em que as linhas vermelha e azul se encontram — são literalmente iguais ao logaritmo escalonado e deslocado. Fora isso, é um bom texto
[010000, 010111]contém 2, 2,25, 2,5, 2,75, 3, 3,25, 3,5 e 3,75Mas as mantissas implícitas pelos logaritmos de base 2 desses números são, respectivamente,
.0000000,.0010101,.0101001,.0111010,.1001010,.1011001,.1100111,.1110100, e, exceto pela primeira, não coincidem com001,010etc. do float. Como os floats no intervalo[2,4)têm espaçamento linear, mas os logaritmos correspondentes não, dá para ver o float, como diz o artigo, como uma aproximação linear por partes do logaritmoEm um gráfico completo haveria
2^32opções dentro do padrão linear por partes, mas o que o texto original desenhou não é esse gráfico completo. Como o texto trata de inteiros de 32 bits e operações com floats IEEE-754 de 32 bits, acho aceitável omitir “discreto” na explicaçãoÉ um bom texto que explica muitos conceitos interessantes, mas a manipulação algébrica em uma seção é surpreendentemente ruim
Na derivação depois de “os passos exatos para ir da primeira forma a esta são muitos, mas incluí todos por completude”, há várias etapas desnecessárias e vários erros de sinal que se anulam entre si. Em especial, ao passar da segunda para a terceira linha, o sinal negativo não é distribuído corretamente. A partir da segunda linha, seria possível começar em
y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2)e chegar muito mais rapidamente ay_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n), com as etapas intermediárias corretas. Para quem entende álgebra, parecem ser apenas passos óbviosO número mágico do famoso trecho de código não é a constante ótima. Usando outra constante, provavelmente dá para reduzir o erro relativo em cerca de 0,5%
Na época talvez fosse difícil encontrar o valor absolutamente ótimo, mas hoje é relativamente fácil. Eu mesmo já caí nessa toca de coelho e tenho um notebook Jupyter para encontrar os números mágicos ótimos de
(1/x^2)e(1/x)O mais interessante para mim neste texto foi o link “How Java's Floating-Point Hurts Everyone Everywhere”: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
O autor é William Kahan, também conhecido como “Old Man of Floating-Point”: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)
JAVAhurte a composição tipográfica é horrível. Parece que usaram algum pacote TeX que aumenta demais o espaçamento entre palavras, e ainda de forma irregular, ou que fizeram OCR de outro documento e entraram espaços extrasHá espaçamentos adicionais estranhos até nas partes em fonte monoespaçada. Foi realmente difícil me concentrar na leitura e, embora eu saiba que não é o caso, quase pareceu um manifesto de um nerd científico
Este vídeo que vi tempos atrás era realmente muito bom: https://www.youtube.com/watch?v=p8u_k2LIZyo