SQRT

O objectivo deste exercício foi criar um algoritmo para encontrar raízes quadradas inteiras usando apenas operações com números inteiros. Desenvolvi inicialmente este método em 1990 para o ZX Spectrum +3 que vinha equipado com um processador Zilog Z80 (compatível com o Intel APX 8080), processador esse que não tem sequer instruções assembly de multiplicação ou divisão!

O desenvolvimento deste documento iniciou-se em 1994, antes de ter fundado a Cynergi, mas não foi terminado até Julho de 2000. Enquanto terminava o documento, encontrei uma melhoria ao algoritmo (a pesquisa binária) que parece ser inovadora e se implementada em hardware, poderia provavelmente gerar uma raiz quadrada de virgula flutuante de 32 bits em 25 ciclos de relógio ou menos, o que é mais rápido do que a ultima geração de processadores Intel (o Itanium, em Agosto de 2000)! Simuladores de FPU (unidade de virgula flutuante, ou floating-point unit, em Inglês) incluídos em sistemas operativos também podem usar este algoritmo para implementar uma raiz quadrada rápida que seja apenas aproximadamente duas vezes mais lenta do que a versão em hardware (num i486).

Dada a sua simplicidade e rapidez, ambos os algoritmos podem também ser usados em

• Sistemas embebidos,
• Máquinas Windows CE,
• Telefones UMTS de 3ª geração,
• Browsers VRML,
• Aceleradores gráficos 3D,
• Calculadoras em relógios,

e quaisquer outros sistemas onde a quantidade de memória é pequena, uma FPU não esteja disponível e/ou onde gráficos 3D (que geralmente necessitam de cálculos de raiz quadrada) sejam necessários. No entanto, e uma vez que o algoritmo da pesquisa binária é mais rápido e barato do que as implementações actuais da raiz quadrada, poderá substitui-las em todas as suas implementações.

Dividi o documento em duas partes, mas deixei-o quase totalmente na mesma ordem em que as ideias foram desenvolvidas. Isto permite que se conheça a sequência de descobertas. Esta Parte I lida com a descoberta empírica de um algoritmo sequencial para encontrar a raiz quadrada, e a sua posterior prova matemática. A Parte II irá lidar com a variante "pesquisa binária" que melhora de forma drástica a velocidade de cálculo da raiz (apesar deste algoritmo sequencial ser o mais pequeno que irá encontrar, se tamanho é o que lhe interessa). Ambos os algoritmos têm a mesma precisão.

O documento completo consiste em:

Parte I

• Descoberta Empírica
• Comparação de Velocidade
• Prova Matemática
• Raiz de Qualquer Índice
• Raiz Cúbica
• Trabalho a Desenvolver

Parte II

• Pesquisa Binária
• Comparação de Velocidade
• Versão em Hardware
• Virgula Flutuante
• Raiz de Qualquer Índice
• Trabalho a Desenvolver
• Conclusão

NOTAS: O endereço "oficial" para este documento é http://www.pedrofreire.com/sqrt. Por favor não crie um atalho directo para o endereço que surge no seu browser — este poderá ser alterado sem aviso prévio. Note ainda que a esta data (Agosto de 2000) não estou ciente da existência ou uso destes algoritmos (pelo menos da vertente "pesquisa binária"), embora não o possa confirmar. Por fim, uso neste documento o ponto para representar a virgula decimal de forma a evitar confusões já que sou obrigado a usá-lo nos trechos de código por definição das linguagens.

Descoberta Empírica

Comecemos pela raiz mais comum (a raiz quadrada ou sqrt(x), do Inglês "square root"). Podemos montar uma tabela de raízes quadradas (onde r é o resultado inteiro):

A partir desta tabela poderá reparar que para números sucessivos, os resultados repetem-se com um padrão, nomeadamente

Com estes valores em mente, podemos agora facilmente criar o seguinte código em C

register unsigned int x, r;

/* … */
x = (x+1) >> 1;    /* init x */
for( r=0; x>r; x-=r++ );

que assume x sem sinal (unsigned). Ao terminar, r tem a raiz inteira, e x tem algum tipo de resto (este e outros restos não serão explorados neste documento).

Experimentemos alguns exemplos para ver se funciona:

O número máximo de iterações para o ciclo for para palavras de n bits é a raiz do maior número que pode ser codificado nesse número de bits, ou seja

`sqrt(2^n - 1) ~~ sqrt(2^n) = 2^(n/2)`

o que, para alguns valores de n, poderá provar que o algoritmo é mais eficiente do que uma conversão inteiro-para-real, seguida da raiz quadrada de virgula flutuante, seguida da conversão real-para-inteiro. Vamos descobrir se tal acontece.

Comparação de Velocidade

Uma optimização em assembly Intel APX/IA do ciclo poderia ser (assumindo AX=x e BX=r, ambos inteiros sem sinal):

	; …
	inc   AX        ;  ciclos i486:    1
	shr   AX, 1     ;                  2
	mov   BX, -1    ;                  1
loop1: 
	inc   BX        ;                  1
	sub   AX, BX    ;                  1
	ja    loop1     ;                3/1

Ao terminar, BX tem o resultado r, e AX algum tipo de resto. Os ponto-e-virgulas foram simplesmente usados para validar os comentários.

O código correspondente usando aritmética de virgula flutuante seria algo como

ctrlword    DW ?

	; Ajustar arredondamento ao inteiro mais proximo
	fstcw  ctrlword
	and    ctrlword, 1111001111111111b
	fldcw  ctrlword

	; …
	fild   WORD PTR [BX]    ;  ciclos i486:    13-16
	fsqrt                   ;                  83-87
	fistp  WORD PTR [BX]    ;                  29-34

assumindo que WORD PTR [BX] tem x. O resultado é devolvido em WORD PTR [BX] (não temos restos em virgula flutuante). Registos de 32 bits podem ser usados em vez dos seus pares de 16 bits em ambos os trechos de código, sem os tornar mais lentos.

Poderá haver uma discrepância entre o resultado devolvido por este código e o devolvido pelo algoritmo. Se um resultado tem uma parte decimal de precisamente .5, o algoritmo irá sempre arredondar o resultado para cima, enquanto a FPU (unidade de virgula flutuante, ou floating-point unit, em Inglês) do i486 irá arrendondá-lo para o número inteiro par mais próximo. Ambos os resultados são correctos, no entanto, e nenhum número inteiro tem uma raiz quadrada com uma parte decimal de .5.

Para comparar a velocidade dos dois métodos, primeiro mostramos que para o algoritmo (chamemo-lhe soft-sqrt), o número de ciclos é de

1+2+1+(1+1+3)(r+1)-2 = 7+5r

onde r é de novo a raiz calculada. Para o método usando a FPU (chamemo-lhe hard-sqrt), temos

Melhor caso:

13+83+29 = 125

Pior caso:

16+87+34 = 137

Logo, a maior gama de raízes que pode ser calculada mais rapidamente pelo soft-sqrt é de

Isto quer dizer que o soft-sqrt é mais rápido do que o hard-sqrt para valores de x entre 0 e 557~676. Cálculos similares para outros processadores APX/IA (x86) geram resultados semelhantes (a melhor gama sendo para o processador i286).

Se precisa de calcular com frequência raízes inteiras na gama acima (que tem n>9, i.e. aproximadamente números de 9 bits), então este algoritmo é-lhe extremamente útil — lembre-se, apenas os valores altos da gama são tão lentos quanto a solução por hardware! Melhor ainda, processadores anteriores ao i486 não têm uma FPU incorporada pelo que a raiz estará provavelmente a ser calculada por um algoritmo lento em software (lento porque tem de levar em conta a virgula flutuante).

A principal desvantagem deste algoritmo é que não tem a precisão da virgula flutuante (devolve o resultado arredondado ao inteiro mais próximo, mas as operações seguintes com este valor podem aumentar o erro).

Prova Matemática

O quadrado de um número r é, naturalmente,

`r^2`

e o quadrado do seu sucessor (número inteiro seguinte) é

`(r+1)^2 = r^2 + 2r + 1`

isto é, está deslocado 2r+1 do r² original. Isto indica-nos a que distância (em termos de x, que é = r²) estará um novo resultado r+1.

Contudo, os valores de r que se alimenta a esta fórmula devem ser escolhidos com cuidado. A fórmula opera com números reais e o meu algoritmo com números inteiros. O meu algoritmo foi desenhado para devolver o resultado da raiz quadrada arredondado para o inteiro mais próximo, o que significa que a alteração de r = 0 para r = 1 por exemplo, irá ocorrer em números reais quando r = 0.5. Isto significa que devemos aplicar 0.5 à fórmula para descobrir quantos resultados de r = 1 existem, já que todos os r = 0.5 a r = 1.49999… serão arredondados para o inteiro 1, e 0.5 é o valor mais baixo neste intervalo. Você quer descobrir a que distância é que r = 0.5 se torna em r = 0.5+1.

Portanto em vez de aplicar

r = 0, 1, 2, 3, 4, 5, …

deveria na realidade aplicar

r = 0, 0.5, 1.5, 2.5, 3.5, 4.5, …

uma vez que 0 é o menor número (maior ou igual a 0) que pode ser arredondado a 0, 0.5 a 1, 1.5 a 2 e assim por diante. A fórmula acima gera assim

o que é obviamente a mesma sucessão de inteiros que foi exibida nos exemplos iniciais. No entanto, se ainda tem dúvidas, então poderá reparar que

como tínhamos assumido.

Q.E.D.

Raiz de Qualquer Índice

Seguindo a mesma linha de raciocínio de acima, podemos ainda resolver raízes de qualquer índice positivo e inteiro (chamemo-lhe i). Para tal, vamos fazer uso do binómio de Newton

`(r+1)^i = sum_(k=0)^i (i!r^(i-k)) / (k!(i-k)!)`

que se reduz a

`(r+1)^2 = (2!r^2)/(0!*2!) + (2!r^1)/(1!*1!) + (2!r^0)/(2!*0!) = r^2 + 2r + 1`

para o caso particular da raiz quadrada (i=2) (e é igual à fórmula a que chegámos acima, verificando ambas). Uma vez que quando o índice do somatório k=0 é

`(i!r^(i-0)) / (0!(i-0)!) = r^i`

o deslocamento (ou distância) pode ser descrito genericamente como

`(r+1)^i - r^i = sum_(k=1)^i (i!r^(i-k)) / (k!(i-k)!)`

e a variação da distância, como

`(((r+1)+1)^i - (r+1)^i) - ((r+1)^i - r^i) = (r+2)^i - 2(r+1)^i + r^i`

ou

`sum_(k=1)^i (i!(r+1)^(i-k)) / (k!(i-k)!) - sum_(k=1)^i (i!r^(i-k)) / (k!(i-k)!) = sum_(k=1)^(i-1) (i!((r+1)^(i-k)-r^(i-k))) / (k!(i-k)!)`

onde o índice superior do somatório é agora i-1 porque quando k=i,

`(i!((r+1)^(i-i)-r^(i-i))) / (i!(i-i)!) = 0`

Se aplicar i=2 a qualquer uma destas fórmulas, elas irão devolver a constante 2, verificando-as e aos resultados já calculados.

Agora deve compreender porque é que usei o binómio de Newton: apesar da primeira destas fórmulas não o mostrar, a segunda demonstra claramente que a variação da distância nunca é de um grau superior a i-2. Porquê i-2 e não i-1? Bom, porque se tomar

`(r+1)^(i-k) - r^(i-k)`

da fórmula com o somatório, verá que esta se assemelha à fórmula do deslocamento ou distância, e a sua resolução mostra que se baixa o resultado por mais um grau.

Como o algoritmo explicado acima necessita da variação da distância para funcionar, altos expoentes (ou graus) implicam um algoritmo lento em altos valores de i. Mas podemos criar código genérico em C para a raiz de qualquer índice:

register int x, s;    /* use double se necessario */
register int r;

/* … */
r = 0;
if( x > 0 )
	{
	s = FD(0.5);
	for( r=1; (x-=s)>0; r++ )
		s += FVD(R-0.5);
	}

Este código só funciona com números positivos: se quiser que funcione com números negativos simplesmente inverta o sinal de r quando x for negativo à entrada, e use abs(x) no interior do código. A operação x-=s anterior à comparação com 0 é a razão de x não ser unsigned. s e r não são unsigned para evitar conversões desnecessárias.

Deve substituir FD(0.5) com o resultado da fórmula do deslocamento ou distância quando r=0.5, e FVD(R-0.5) pela fórmula da variação da distância, substituindo r por r-0.5. Deveria também tornar as variáveis x e s double (reais) em vez de int (inteiras) se quaisquer destas substituições devolve números não-inteiros. Isto irá desacelerar o código, mas vai ter de o optimizar à mão, de qualquer forma. x deverá ter o valor para o qual se quer calcular a raiz quadrada, ao entrar no código, e r irá ter a raiz calculada à saída.

Por exemplo, para a raiz quadrada, deveria fazer FD(0.5) = 2 e FVD(R-0.5) = 2 uma vez que é constante. x e s mantêm-se int.

Qual é a rapidez deste código, em comparação com a correspondente operação de hardware? Bem, isso depende de quatro factores:

1. Com que rapidez consegue o seu CPU calcular essa raiz em virgula flutuante;
2. Precisa de doubles em vez de ints;
3. Qual a rapidez da multiplicação inteira no seu CPU;
4. Qual a qualidade do código optimizado que consegue gerar.

O primeiro factor diz-lhe se todo este trabalho vale de facto a pena. Na maioria dos CPUs, cálculos de virgula flutuante levam dezenas (se não centenas) de ciclos de relógio a completar, e a única raiz que fazem é a quadrada (i=2). Outras raízes necessitam de ser calculadas usando logaritmos o que torna essas operações ainda mais lentas e este algoritmo ainda mais rápido (comparativamente). Por outro lado, para valores altos de i, o código do ciclo interno deste algoritmo abrange mais números (em termos de x — isto quer dizer que s cresce mais rapidamente), o que quer dizer que o código chega mais rapidamente a uma solução. Mas se o seu CPU consegue executar uma raiz de forma muito rápida, todo o trabalho aqui descrito é inútil — use-a!

O segundo factor (usar double em vez de int) pode ser evitado se souber suficiente matemática e souber como optimizar código. Por exemplo, um valor de 3.25 para FD(0.5) pode ser transformado em 3, e um novo contador inicializado para subtrair 1 de x a cada 4 iterações. Dependendo da rapidez de execução de adições (e subtracções) em virgula flutuante (reais) no seu CPU, isto poderá ou não valer a pena.

O terceiro factor determina até que ponto podemos subir em termos de i antes do algoritmo deixar de ser eficiente. Alguns CPUs executam uma multiplicação inteira em um ciclo de relógio. Para esses, de facto não vai importar muito o grau de FVD(R-0.5) já que apenas precisará de fazer uma sequência de multiplicações, e essas são baratas o suficiente (em termos de velocidade) nesses CPUs.

Finalmente, o ultimo factor está inteiramente dependente de si. Para i=2, por exemplo, o código acima torna-se no código que exibi ao início deste documento — sem ifs e com um ciclo muito apertado. Só a experiência o pode ajudar a este ponto. — não é o objectivo deste documento dar-lhe uma série de algoritmos optimizados para vários valores de i, especialmente uma vez que tal optimização seria feita em assembly e variantes dessa linguagem abundam. Irei no entanto focar-nos na raiz cúbica para o preparar para tal análise.

Raiz Cúbica

Como exercício e exemplo, vamos criar um algoritmo para calcular a raiz cúbica (i=3). Vamos aplicar o valor de i à fórmula do deslocamento ou distância:

`sum_(k=1)^3 (3!r^(3-k)) / (k!(3-k)!) = (3!r^2)/(1!*2!) + (3!r^1)/(2!*1!) + (3!r^0)/(3!*0!) = 3r^2 + 3r + 1`

Quando faz r=0.5, obtém o valor 3.25 desta fórmula. Isto quer dizer que FD(0.5) = 3.25. Agora fazemos o mesmo à fórmula da variação da distância:

`sum_(k=1)^2 (3!((r+1)^(3-k)-r^(3-k))) / (k!(3-k)!) = (3!((r+1)^2-r^2))/(1!*2!) + (3!((r+1)^1-r^1))/(2!*1!) = 6r + 6`

Aqui, quando faz r=r-0.5, obtém 6(r-0.5)+6 = 6r+3. Isto quer dizer que FVD(R-0.5) = 6r+3. Isto dá-nos o algoritmo genérico

register double x, s;    /* double *e'* necessario */
register int r;

/* … */
r = 0;
if( x > 0 )
	{
	s = 3.25;
	for( r=1; (x-=s)>0; r++ )
		s += 6*r+3;
	}

Mais uma vez, este código só funciona com números positivos: se quiser que funcione com números negativos simplesmente inverta o sinal de r quando x for negativo à entrada, e use abs(x) no interior do código. Tivemos de usar x e s como doubles porque apesar de FVD(R-0.5) nunca devolver um número não-inteiro, o valor inicial FD(0.5) não é inteiro. Portanto, o nosso primeiro passo em optimizar este código é livramo-nos dos valores em double, tornando todo o ciclo inteiro.

Note que os passos seguintes só se aplicam a esta raiz em particular, e podem não se aplicar a outras. Para aqueles não familiarizados com a linguagem C, note-se ainda que x-=s significa x=x-s, o que é uma expressão válida devolvendo o resultado da subtracção.

A parte decimal de .25 de s nunca se altera: s é incrementado por 6r+3, o que sempre devolve inteiros se alimentado com valores inteiros de r (o que é). A parte decimal presente no valor inicial de s tem no entanto um impacto em cada iteração de x, no entanto.

Uma vez que .25 = ¼, torna-se 1 em cada 4 iterações do ciclo. Isto significa que se deveria subtrair 1 de x quando a sua parte decimal for .00, antes da subtracção de s. Porquê antes da subtracção? Bem, porque precisaríamos de verificar se este decremento faz x<0, e isto nunca acontece nesse ponto, poupando-nos outra comparação.

A parte decimal de x é .00 ao início da primeira iteração do ciclo e ao início de cada 4 iterações após essa. Vamos ver um exemplo que começa em x=20 e s=3.25.

Portanto você decrementa x antes da primeira iteração do ciclo, e volta a fazê-lo antes de x-=s quando os dois bits menos significativos de r são iguais a um certo valor. Isto elimina a necessidade de um contador — afinal de contas, r comporta-se como um contador, e os dois bits menos significativos de um contador repetem-se a cada 4 contagens, automaticamente. Mas que valor escolher?

Bem, o valor na coluna iteração na tabela acima é o mesmo que r dentro do ciclo (r é incrementado ao final do ciclo, após tudo o resto). Aqui poderá ver que precisa de decrementar x antes da iteração 1 e antes da iteração 5. Isto significa iterações 0 e 4, ou quando os 2 bits menos significativos de r=00b.

Isto ajuda a manter sincronizados os valores de s e x. Comparar ambos (na subtracção x-=s) é agora trivial excepto quando são iguais. Ao contrário do código com double, não podemos agora assumir que x-=s igual a 0 seja um sinal para terminar o ciclo, porque mesmo que as suas partes inteiras sejam iguais, as suas partes decimais podem não o ser, pelo que precisamos de uma forma de descobrir qual era a parte decimal de x logo após a subtracção. Para tal, vamos analisar melhor a iteração do ciclo nesse ponto:

Isto é válido para as primeiras 4 iterações do ciclo, e para cada 4 após essas. Tal como esta tabela mostra, só quando os 2 bits menos significativos de r são 00b é que o ciclo terminaria na realidade, caso contrário ocorreria mais uma iteração. Portanto tudo o que o código tem de fazer nestas situações (mais uma iteração) é incrementar r antes de terminar.

Vamos ver uma nova versão do código que leva tudo isto em atenção:

register int x, s;
register int r;

/* … */
r = 0;
if( x > 0 )
	{
	x--;
	s = 3;
	for( r=1; (x-=s)>=0; r++ )
		{
		if( x==0 )
			{
			if( (r & 3) != 0 )
				r++;
			break;
			}
		s += 6*r + 3;
		if( (r & 3) == 0 )    /* 3 == 11b */
			x--;
		}
	}

Agora todas as operações são do tipo int, mas o código ainda só funciona com números positivos. Existe uma comparação dupla de x com 0, mas isso pode ser facilmente corrigido em assembly. O que podemos ainda melhorar é a multiplicação: esta pode ser removida, tornando o código mais rápido. Isto pode ser feito com um contador que se inicializa a 6 e é incrementado por 6 a cada iteração. Irá então conter 6r já que r é um contador que se inicializa a 1 e é incrementado por 1 a cada iteração. Isto substitui uma multiplicação por uma adição, uma mudança muito bem-vinda.

Poderia também argumentar-se que o contador para r já não seria necessário uma vez que podemos determinar r dividindo o contador da multiplicação por 6 ao final do ciclo. Isto é verdade e ainda poderia testar os bits 1 e 2 deste contador por 00b em vez de testar os bits 0 e 1 r pelo mesmo valor (não irei demonstrar isto). Mas a operação de divisão é muito cara (24 a 40 ciclos de relógio) versus aquilo que poupa (1 ciclo de relógio por iteração), se quiser usar esta "soft-raiz-cubica" em vez de uma "hard-raiz-cubica". Em alguns CPUs (e de forma mais geral, se não tiver uma FPU em hardware), isto poderá compensar, no entanto.

Uma vez que o código C para a ultima alteração é trivial, vamos passar directamente a uma versão deste algoritmo em assembly (assumindo AX=x e BX=r, ambos inteiros sem sinal):

	; …
	mov   BX, 0      ;  BX = r         ciclos i486:        1
	sub   AX, 1      ;                                     1
	jb    done       ;                                   3/1
	mov   CX, 3      ;  CX = s                             1
	mov   DX, 0      ;  DX = contador da multiplicacao     1
loop1:
	inc   BX         ;                                     1
	sub   AX, CX     ;                                     1
	jb    done       ;                                   3/1
	je    equal      ;                                   3/1
	add   DX, 6      ;                                     1
	add   CX, DX     ;                                     1
	add   CX, 3      ;                                     1
	test  BX, 11b    ;                                     1
	jnz   loop1      ;                                   3/1
	dec   AX         ;                                     1
	jmp   loop1      ;                                     3
equal:
	test  BX, 11b    ;                                     1
	jz    done       ;                                   3/1
	inc   BX         ;                                     1
done:

À saída, BX tem o resultado, r. Este código está optimizado para velocidade num i486, e não para tamanho. Por outro lado, as operações test não estão disponíveis num processador anterior ao i386, o que quereria dizer que teria de encontrar uma alternativa às mesmas usando uma cópia temporária de BX e um and.

Agora vamos escrever o código de virgula flutuante para o mesmo processador. A família de processadores Intel APX/IA-16/32 (i.e., o x86 até ao Pentium-III) não tem uma operação hardware para a raiz de qualquer índice. Poderíamos então usar

`root(i) x = x^(1/i)`

mas eles não têm sequer uma operação para elevar qualquer base a qualquer expoente! Eles têm no entanto logaritmos selectos e bases selectas elevadas a expoentes pelo que pode usar a regra

`x = b^(log_(b) x)`

para construir a fórmula

`x^(1/i) = b^(log_(b) x^(1/i))`

`x^(1/i) = b^(1/i log_(b) x)`

que devolve a raiz de índice i de x, dada qualquer base real, b.

Computadores trabalham na base 2, pelo que faz sentido que logaritmos e expoentes que eles calculam sejam na base 2. Isto é verdade para o i486 (pelo que fazemos b=2). O i486 tem as seguintes operações relevantes:

Cada linha é uma operação do CPU: esta operação executa mais do que uma operação matemática de cada vez, o que pode ou não ajudar-nos. Algumas são limitadas na gama de números aceites (as variáveis que não são exibidas nas limitações, não são limitadas). Mas algo é óbvio — não existe operação para elevar 2 (ou outra base qualquer) a um expoente arbitrário, e nós precisamos disso! Felizmente você pode usar

`b^m * b^n = b^(m+n)`

para resolver esse problema usando as duas primeiras operações exibidas acima. Para tal faça m=int(z) (sendo z o expoente a que quer elevar 2), e n=z-int(z). Você atribui estes valores aos expoentes de cada uma das primeiras duas operações exibidas acima, respectivamente, e é tudo! Uma vez que int(z) é sempre um inteiro e z-int(z) está sempre entre -1 e 1 (exclusivamente), independentemente do método de arredondamento, este é o passo final para traduzir a nossa raiz para:

`z = 1/i * log_(2) x`

`r = ((2^(z-"int"(z))-1)+1) * 2^("int"(z))`

Vamos usar o arredondamento para o inteiro mais próximo no código, já que precisamos desse tipo de arredondamento para o resultado final, e não queremos perder tempo a trocar de arredondamento a meio de um cálculo. Isto quer dizer que z-int(z) estará na realidade entre -0.5 e 0.5 (inclusivamente). O código final é então:

ctrlword    DW ?
inv_i       REAL8 1/i    ; substitua i ao "assemblar"/montar
one         REAL8 1

	; Ajustar arredondamento ao inteiro mais proximo
	fstcw  ctrlword
	and    ctrlword, 1111001111111111b
	fldcw  ctrlword

	; …
	fld    inv_i            ;                   ciclos i486:       3
	fild   WORD PTR [BX]    ;                                 13- 16
	fyl2x                   ;  ST(1) = (1/i) * log2 [BX]     196-329
	fld    ST               ;  dup ST: ST(1) now = ST              4
	frndint                 ;  ST = int(ST)                   21- 30
	fxch                    ;  ST <-> ST(1)                        4
	fsub   ST, ST(1)        ;  ST = ST - ST(1)                 8- 20 (*)
	f2xm1                   ;  ST = 2^ST - 1                 140-279
	fadd   one              ;  ST++                            8- 20
	fscale                  ;  ST = ST * 2^ST(1)              30- 32
	fistp  WORD PTR [BX]    ;                                 29- 34
	ffree  ST               ;  (veja a proxima linha)              3
	fincstp                 ;  "pop" do velho ST(1)                3
	;
	;                 (*) neste ponto (apos o fsub):
	;                     se z = (1/i) * log2 [BX], entao
	;                        ST(1) = int(z)
	;                        ST    = z - int(z)
	;                        ST esta' entre [-0.5 .. 0.5]

A operação f2xm1 só opera no intervalo [0 .. 0.5] numa FPU 8087. Neste processador você precisa de verificar se ST era negativo. Se sim, então precisa de alimentar o seu valor absoluto a f2xm1 e encontrar o inverso do resultado. As ultimas duas operações são necessárias para remover da pilha um valor temporário que não é removido pela ultima operação, fscale. Se ignora os avisos que a FPU lhe pode dar sobre a pilha estar cheia, pode ignorar as operações e poupar algum tempo.

A parte "divertida" deste código é que na realidade eleva o valor inteiro guardado em WORD PTR [BX] ao expoente real guardado em inv_i. Se guardasse nesta ultima localização o valor 8.25, por exemplo, este ainda seria o código optimizado para elevar um inteiro em WORD PTR [BX] a esse expoente. Para calcular uma raiz de índice i, terá de guardar o inverso de i (tal como exibido) nessa localização, antes da primeira execução do código. Isto pode ser feito pelo próprio assembler. Isto também significa que este código se aplica a toda as outras raízes que não a raiz quadrada.

Vamos agora comparar a velocidade do algoritmo versus este código baseado em hardware.

Começamos por mostrar que para a soft-raiz-cubica, o número de ciclos de relógio é de 5 se r=0, ou

5×1+(8×1+3)(r-1)+(-2+1+3)(¼)(r-1)+5 = 11.5r - 1.5

se r>0, onde r é de novo a raiz calculada. O número de ciclos de relógio das ultimas 3 operações (entre as etiquetas — labels — "igual" e "feito"), que ocorre a cada ¼ (25%) dos possíveis valores finais de r, ou cada

`(1/4 * sqrt(2^n))/(2^n) = 2^(-n/2-2)`

cálculos aleatórios da raiz cúbica para números com n bits, é negligenciável e portanto não adicionado a esta fórmula (e.g.: para n=16 você tem essas operações a serem executadas menos de 0.1% vezes, e o valor a ser adicionado a esta fórmula seria inferior a 0.004). Para a hard-raiz-cubica, você tem

Melhor caso:

3+13+196+…+29+3+3 = 462

Pior caso:

3+16+329+…+34+3+3 = 777

Assim, a maior gama de raízes que pode ser calculada mais rapidamente pela soft-raiz-cubica é de

Isto significa que a soft-raiz-cubica é mais rápida do que a hard-raiz-cubica para valores de x entre 0 e 1624~4583. Cálculos similares para outros processadores APX/IA (x86) geram resultados semelhantes (a melhor gama sendo para o processador i286).

Se precisa de calcular com frequência raízes inteiras na gama acima (que tem n>10, i.e. aproximadamente números de 10 a 13 bits), então este algoritmo interessa-lhe. Ele é muito mais apelativo do que o soft-sqrt era para raízes quadradas!

Trabalho a Desenvolver

Este trabalho ainda pode ser estendido em pelo menos duas direcções: precisão e velocidade.

Em termos de precisão, você pode fazer os cálculos necessários para de alguma forma fazer uso dos valores devolvidos (para além de r) como restos, de forma a fazer cálculos mais precisos. A aplicação disto é limitada, no entanto.

Também pode "afinar" o algoritmo de forma a que quando a parte decimal do resultado real r for .5, ele arredonde o resultado ao número inteiro par mais próximo e não sempre para cima como faz agora. Tal como já foi dito, raízes quadradas de números inteiros nunca têm partes decimais de precisamente .5, pelo que isto não seria mais do que um exercício para soft-sqrt. Para outras raízes, tal precisão é na minha opinião inútil já que estamos a falar de números inteiros afinal de contas, e a precisão já foi pela janela fora, mas tornaria o algoritmo talvez um pouco mais preciso.

Em termos de velocidade, poderia sempre tentar usar uma pesquisa binária (o tipo de pesquisa que usa quando está à procura de uma palavra num dicionário), se esta se provar rápida o suficiente para encontrar raízes, ou simplesmente inverter a direcção de pesquisa pela raiz (começando pelos valores mais altos de x e r e pesquisando para trás). Isto significaria que as raízes maiores são pesquisadas primeiro, pelo que uma maior gama de valores de x são cobertos nas primeiras iterações, e em média uma raiz é encontrada mais cedo. O pior caso ainda é uma pesquisa total, no entanto.

Por exemplo, uma optimização em assembly Intel APX/IA desta variante poderia ser (assumindo AX=x e BX=r, ambos inteiros sem sinal, e n=16, onde n significa o número de bits a serem calculados):

	; …
	not   AX         ;  AX = (2^n)-1-AX       ciclos i486:   1
	mov   BX, 256    ;  BX = 2^(n/2)                         1
	sub   AX, 254    ;  AX = AX - ( (2^(n/2))-2 )            1
	jbe   done       ;                                     3/1
	inc   AX         ;                                       1
	shr   AX, 1      ;                                       2
loop1:
	dec   BX         ;                                       1
	sub   AX, BX     ;                                       1
	ja    loop1      ;                                     3/1
done:

Para n=16, este algoritmo tem um número de ciclos de relógio de 6 se r=256, ou

7+5×(256-r)-2 = 1285-5r

se r<256, o que significa que consegue calcular raízes quadradas cujos resultados sejam de 230~232 a 256 durante o mesmo tempo que leva a calcular hard-sqrt. Isto é uma gama de 12636 valores de x (19.3%), comparado com apenas 677 (1%) para o algoritmo já exibido. O pior caso é idêntico em ambos, no entanto (ou algo pior neste algoritmo, mas a diferença é fixa e não proporcional a r).

Iremos ainda prestar mais atenção à variante de pesquisa binária na Parte II deste documento, mas outro trabalho é deixado para o leitor.