Introdução
Algoritmos derivados de 0x5F3759DF
começam usando os ditos "números mágicos" para em forma de inteiro formar o código binário do float
que aproximará a solução, podendo depois refinar a precisão com operações de ponto flutuante, normalmente utilizado o método de newton até três vezes (na maioria das vezes, só uma ou duas vezes).
Além disso, para ter melhor desempenho se deixa de tratar casos específicos (argumentos próximos de zero, negativos, nan, etc) e foca-se em casos mais comuns, reservando esses tratamentos à aplicação de acordo com as necessidades.
Algoritmos sem aproximações refinadas
A forma mais rápida de calcular sqrt
e invsqrt
somente faz uma soma/subtração e um deslocamento de bits com o código binário do argumento em forma de inteiro (eu prefiro sem sinal), sendo assim muito rápido, porém pouco preciso.
O algoritmo a seguir calcula sqrt
de um radicando (matematicamente, radicando é o argumento de uma raiz quadrada, cúbica, qualquer raiz em qualquer índice) com erro relativo abaixo de 3.475% em condições normais. Quando o radicando é zero, o resultado é 7.933e-20, resultando nesse mesmo erro absoluto neste caso. Por conveniência, vou assumir que a execução desse algoritmo leva um sqrtTime
.
uint code=( 0x3F769E5C + *(uint*)&radicand )>>1 ;
return *(float*)&code ;
Já o seguinte calcula invsqrt
de um radicando com erro relativo menor que 3.422%. Quando o radicando é infinito, o resultado é 5.239e-20, resultando nesse mesmo erro absoluto neste caso. Custa cerca de um sqrtTime
mesmo.
uint code=( 0xBE6EC85F - *(uint*)&radicand )>>1 ;
return *(float*)&code ;
Algoritmos com uma iteração de método de newton
De fato, erros relativos acima de 1% são inconvenientes. Devido a isso, ao menos uma vez se aplica sobre esses valores uma iteração do método de newton adaptado à função específica.
No caso da raiz quadrada, isso consiste em calcular NewApproach=( OldApproach + Radicand/OldApproach )/2
, já no caso da inversa, NewApproach = OldApproach*( 3 - OldApproach*OldApproach*Radicand )/2
.
Além disso, já que [( √r1=Sqrt )]&&[( r2^(-1/2)=Sqrt )]=>[( r1=1/r2 )]
, também pode-se aplicar na inversa a fórmula NewApproach=( OldApproach + 1/( Radicand*OldApproach ) )/2
, que é uma suave modificação do método de newton para refinar a raiz quadrada.
Sendo assim, temos a seguir um algoritmo de sqrt
que tem erro relativo menor que 6.011e-4 e, segundo meus testes, custa pelo menos 2.8*sqrtTime
.
uint code=( 0x3F76CF5E + *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
return 0.5f*( root+radicand/root ) ;
Agora um algoritmo de invsqrt
que tem erro relativo menor que 1.752e-3 e custa pelo menos 2.5*sqrtTime
.
uint code=( 0xBE6EB50D - *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
return root*( 1.5f + root*root*( radicand*-0.5f ) ) ;
E finalmente um algoritmo de invsqrt
que tem erro relativo menor que 5.895e-4 e custa pelo menos 3.3*sqrtTime
.
uint code=( 0xBE6EB50D - *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
return 0.5f*( root+1/( radicand*root ) ) ;
Algoritmos com duas iterações de método de newton
Normalmente, considera-se que mais de duas iterações têm custo/benefício alto, por isso vamos parar por aqui com as iterações. Com duas delas, percebe-se que as vezes é possível simplificar operações e também que invsqrt
pode utilizar duas vezes uma das duas fórmulas ou utilizar uma vez cada. Se utilizar uma de cada, em questão de desempenho pouca diferença faz qual vem primeiro, mas a precisão é maior quando primeiramente se utiliza a mais precisa.
Temos a seguir um algoritmo de sqrt
que tem erro relativo menor que 1.805e-7 e, segundo meus testes, custa pelo menos 3.9*sqrtTime
.
uint code = ( 0x3F76CF5E + *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
root += radicand/root ;
return 0.25f*root + radicand/root ;
Agora temos a seguir um algoritmo de invsqrt
que tem erro relativo menor que 4.598e-6 e, segundo meus testes, custa pelo menos 4.3*sqrtTime
. Não se sabe a razão desse algoritmo ser mais lento que o anterior e, ainda assim, a retirada de uma iteração de cada um o torna mais veloz. Pode ser consequência de simplificações no outro.
uint code = ( 0xBE6EB50D - *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
root *= 1.5f + root*root*( radicand*=-0.5f ) ;
return root*( 1.5f + root*root*radicand ) ;
Agora, outro algoritmo de invsqrt
. Este tem erro relativo menor que 5.213e-7 e, segundo meus testes, custa pelo menos 4.4*sqrtTime
.
uint code = ( 0xBE6F02E3 - *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
root = root+1/( radicand*root ) ;
return root*( 0.75f + root*root*( radicand*-0.0625f ) ) ;
E finalmente o último algoritmo de invsqrt
com iterações. Este tem erro relativo menor que 1.737e-7 e, segundo meus testes, custa pelo menos 5.2*sqrtTime
. Percebe-se em relação ao anterior uma estranhamente enorme perda de desempenho e ínfimo ganho de precisão, o que torna este algoritmo menos recomendado.
uint code = ( 0xBE6F02E3 - *(uint*)&radicand )>>1 ;
float root = *(float*)&code ;
root = root+1/( radicand*root ) ;
return 0.25f*( root+4/( radicand*root ) ) ;
Algoritmos com tabelas de números mágicos
Para finalizar, consideremos a possibilidade de refinar a precisão sem operações com float
mas sim com uma seleção de distintas fórmulas na determinação do código binário via tabela de números mágicos em array. Pode-se fazer tabelas de diversos tamanhos, mas considera-se adequado o uso de tabelas com 1024 números inteiros sem sinal de 32 bits.
Os dois algoritmos a seguir que desenvolvi pegam uma parte do código binário do radicando e utiliza como índice de tabelas de termos próximos de 0x3F769E5C
e 0xBE6EC85F
(utilizados nos dois primeiros algoritmos, um de sqrt
e outro de invsqrt
) e também como índice de uma tabela de fatores que ajustam a derivada da função.
Repare que os fatores são usados num produto do tipo Factor*CodePiece/2^32
que envolvem números inteiros de 32 bits, o que em bons compiladores é facilmente simplificado em proveito das instruções de máquina.
Primeiramente, o algoritmo de sqrt
com erro abaixo de 2.008e-7 e custo de pelo menos 2.5*sqrtTime
segundo as medidas que fiz. Faz uso das tabelas SqrtCodeTerms
e SqrtCodeFactors
.
uint code=*(uint*)&radicand , piece=code&0x00FFFFFF , index=piece>>14 ;
code=( (uint)(( (int)SqrtCodeFactors[index]*(long)piece )>>32)+SqrtCodeTerms[index]+code )>>1 ;
return *(float*)&code ;
E o algoritmo de invsqrt
com erro abaixo de 4.097e-7 e custo de pelo menos 2.7*sqrtTime
. Não se sabe a razão da diferença de desempenho com o algoritmo anterior. Faz uso das tabelas InvSqrtCodeTerms
e InvSqrtCodeFactors
.
uint code=*(uint*)&radicand , piece=code&0x00FFFFFF , index=piece>>14 ;
code=(uint)(( (int)InvSqrtCodeFactors[index]*(long)piece )>>32)+(( InvSqrtCodeTerms[index]-code )>>1) ;
return *(float*)&code ;
Para finalizar, essas tabelas foram obtidas por meio de um programa matemático da maplesoft. Com o maple, programou-se o código que calculava para cada índice os valores mais adequados para os algoritmos e esses valores foram impressos. As tabelas estão nos links a seguir.
Sqrt: SqrtCodeTerms e SqrtCodeFactors
InvSqrt: InvSqrtCodeTerms e InvSqrtCodeFactors