SQRT

I&D: Raiz Quadrada em Assembly | Parte II

Esta é a segunda parte do algoritmo de raiz quadrada. Foi desenvolvida em um mês durante as fases finais de terminar este site e rever este documento para publicação.

Enquanto que a Parte I lidava com a descoberta empírica de um algoritmo sequencial para encontrar a raiz quadrada, e a sua posterior prova matemática, esta Parte II lida com a variante "pesquisa binária" que melhora de forma drástica a velocidade de cálculo da raiz (apesar do 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

Parte II

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.

Pesquisa Binária

Uma pesquisa binária funciona de forma semelhante à pesquisa por uma palavra num dicionário. Quando procura por "Verde", você não abre o dicionário na primeira página e começa a procurar a sua palavra por todo o dicionário. Em vez disso você tira partido do conhecimento de que as palavras estão ordenadas alfabeticamente.

Isto é como descreveria a pesquisa a outra pessoa: Primeiro abra o dicionário em metade, e olhe para a primeira palavra que encontra. Compare-a com a palavra por que está a procurar. Se for igual, encontrou-a. Se não, divida uma das duas metades do dicionário de novo em metade (a metade direita se a palavra por que estava a procurar está depois da que encontrou, a metade esquerda se antes, alfabeticamente), e repita o processo.

Lembra-se do jogo das advinhas, da infância? "Advinhem um número". "Estou a pensar num número entre 1 e 100 - advinhem qual é em 10 tentativas". "Só vos posso dizer se o número que advinham está acima ou abaixo do número em que pensei".

A melhor estratégia para ganhar neste jogo é exactamente a mesma. Você começa por advinhar 50, e depois passa a metade disso, depois metade, e assim por diante, baseado nas respostas que recebe. Se fizer isto irá sempre encontrar o número em menos de 8 tentativas! Isto é chamada uma pesquisa binária.

Para o algoritmo da raiz quadrada podemos fazer uma pesquisa semelhante (porque quer a raiz quadrada, quer o quadrado são operações contínuas e crescentes, logo "ordenadas"): primeiro começamos em metade da raiz máxima que pode ser calculada. Para números de n bits, isto é

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

Por exemplo, para n=16, metade da raiz máxima é 128. Dependendo de como o valor original a encontrar a raiz, x, se compara com o quadrado de 128, você corta uma das duas metades de novo em metade. Um-quarto da raiz máxima (ou metade de 128) é 64, pelo que se deveria escolher a metade inferior, você compara x com o quadrado de 64, se deveria escolher a metade superior, você compara x com o quadrado de 192.

Como 64 está precisamente no meio da metade inferior, pode ser obtido por 0+64 ou 128-64, e como 192 está precisamente no meio da metade superior, pode ser obtido por 128+64 ou 256-64 (i.e., você obtém estes valores ao adicionar um-quarto da raiz máxima à borda inferior da metade, ou ao subtrair esse valor da borda superior, como pode ser visto no gráfico seguinte).

[0------64-----128----192----256]; 192->256 = 128->192 = etc. = 64

Quando você quer dividir o próximo quarto em metade, irá precisar de usar um-oitavo da raiz máxima (32). No próximo passo, um-dezasseis-avos (16) e assim por diante. A forma mais simples de obter o valor (metade da secção para onde se quer mover) a partir destes valores parciais da raiz máxima, é adicioná-los ao valor actual (se quiser a secção superior) ou subtraí-los (se quiser a secção inferior). Por exemplo, se estivesse em 192 e x fosse agora inferior ao quadrado de 192, então o valor seguinte com que quer comparar x é com o quadrado de 192-32, ou 160.

Tudo isto está muito bem, e o quadrado do primeiro valor é relativamente fácil de calcular — afinal de contas,

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

e para cada das metades seguintes (por exemplo, se a secção escolhida fosse sempre a inferior),

`(2^(n/2-2))^2 = 2^(n-4)`

`(2^(n/2-3))^2 = 2^(n-6)`

`(2^(n/2-4))^2 = 2^(n-8)`

…

`(2^(n/2-k))^2 = 2^(n-2k)`

o que significa que o quadrado destas potências de 2 pode ser calculado consecutivamente fazendo uma deslocação (shift) do primeiro quadrado duas vezes para a direita, tal como as raízes quadradas respectivas (as metades da raiz máxima) podem ser calculadas fazendo uma deslocação da primeira metade uma vez para a direita. Mas e se precisarmos de seleccionar uma secção superior? Já lá vamos.

Até agora só temos deslocações, adições, subtracções e comparações inteiras, o que é óptimo para os nossos propósitos. Mas nem todo o método está pronto. Temos uma forma de aproximar as raízes, mas não de aproximar os quadrados. Por exemplo, se altera 128 para 192 adicionando 64, o que deve fazer com os respectivos quadrados? Será que se deve calcular o quadrado fazendo uma multiplicação?

Bem, não necessariamente (o que é bom, dado que isso seria lento!). É que, se a é o actual valor em pesquisa (192 por exemplo), e b a metade da raiz máxima que quer adicionar ou subtrair dela (32 neste exemplo), então o que quer encontrar é

`(a+-b)^2`

e você tem

`(a+-b)^2 = a^2 +- 2ab + b^2`

o que quer dizer que o deslocamento ou distância

`+-2ab + b^2`

deveria ser adicionado ao quadrado actual, sendo o sinal de 2ab idêntico à operação que se faz com as raízes (adição / subtracção).

Aqui temos 3 multiplicações, mas vamos vê-las mais de perto: a primeira é uma multiplicação por 2, o que pode ser substituído por uma operação de deslocação. A seguir multiplica-se por b, que é uma potência de 2, portanto também esta multiplicação pode ser substituída por uma operação de deslocação. Finalmente, b² é simples de calcular exactamente porque b é uma potência de 2, tal como já vimos!

Por outras palavras, o deslocamento ou distância é na realidade

`a * 2^(n/2-k+1) + 2^(n-2k)`

o que, para um n par e k>0 inteiro, devolve sempre expoentes positivos e inteiros de 2, o que quer dizer que todas as operações são deslocações rápidas e adições, e não multiplicações! Isto também prova matematicamente os passos usados no método.

Isto permite-nos criar a primeira versão em C do código deste ciclo: 

#include <math.h>
#include <limits.h>
/* … */
register unsigned long int x, r, h, sq, x2;
register int i;    /* nao vale a pena serem grandes como x, etc */

/* … */
sq = (ULONG_MAX >> 2) + 1;    /* (2^(n/2-1))^2 = 2^(n-2) */
r  = sqrt(sq);                /* 2^(n/2-1) */
i  = log(r) / log(2);         /* n/2-1 */
/* nota: recodifique as 2 linhas acima de forma a que os calculos
   sejam feitos ou por voce^ ou pelo compilador durante a compilacao.
   compiladores inteligentes *podem* ja' fazer isso */
h  = r;
x2 = sq;
do	{
	sq >>= 2;        /* sq = b^2 */
	h >>= 1;         /* h = b */
	if( x >= x2 )    /* x == x2 pode tomar qualquer ramo */
		{
		x2 = x2 + (r << i) + sq;    /* (r << i) == 2ab */
		r += h;
		}
	else
		{
		x2 = x2 - (r << i) + sq;    /* (r << i) == 2ab */
		r -= h;
		}
	}
	while( --i != 0 );
/* nota: x esta' entre x2 e um x2 previamente testado. uma vez
   que na ultima iteracao h=1, ou r e' a raiz que queremos, ou
   um dos seus valores adjacentes o e' */
if( x <= (x2-r) )
	r--;
else if( x > (x2+r) )
	r++;

Aqui, a variável i tem uma função similar ao k da fórmula acima, apesar de não serem idênticas em termos de valores que tomam. As ultimas 4 linhas do código são um produto da análise feita na Parte I deste documento. Ao entrar no código, x deveria ter o valor a encontrar a raiz quadrada, e ao sair, r irá ter o resultado inteiro, arredondado ao inteiro mais próximo.

O que poderá verificar de imediato é que ao contrário dos antigos algoritmos sequenciais que levavam

`sqrt(x)-1` a `sqrt(x)`

iterações a terminar, o que significa que o número de iterações tinha uma gama de

`0 <= sqrt(x)-1 <= 2^(n/2)-1` a `0 <= sqrt(x) <= 2^(n/2)`

a pesquisa binária tem um número fixo de iterações de

`(n/2)-1`

Isto é uma grande melhoria sobre os algoritmos anteriores (embora estes sejam ainda mais rápidos para uma muito pequena gama de raízes). Por exemplo, mesmo para um valor baixo de n (por exemplo n=16), os outros algoritmos levavam entre 0 e 256 iterações para encontrar a raiz, enquanto este leva exactamente 7 iterações para encontrar essa mesma raiz!

Pela mesma razão não se usa x==x2 como condição para sair do ciclo já que só melhora a velocidade para

`2^(n/2-1) - 1`

valores de x e piora para todos os outros!

O próximo passo na optimização pede pela sua atenção nas deslocações. Vê como a única utilização de r no ciclo é com uma deslocação para a esquerda por um número de dígitos que decresce a cada iteração, e é actualizado por um valor (h) que leva uma deslocação para a direita por 1 a cada iteração? Se fosse começar o ciclo com r totalmente deslocado para a esquerda, e se lhe fosse aplicada uma deslocação para a direita por 1 a cada iteração, h teria de ser deslocado para a direita por 2 (a mesma direcção e valor da deslocação de sq). Vê como h e i não são usados em mais nenhum outro lado no código?

Isto permite-nos criar o seguinte código optimizado em C para este ciclo, onde h já não tem o significado de "metade da raiz" ("half"), mas antes de valor "auxiliar" ("helper") que junta as funções do antigo h e sq. 

#include <limits.h>
/* … */
register unsigned long int x, r, h, x2;

/* … */
r  = (ULONG_MAX >> 2) + 1;    /* (2^(n/2-1)) << (n/2-1) */
h  = (ULONG_MAX >> 3) + 1;    /* (2^(n/2-1)) << (n/2-2) */
x2 = (ULONG_MAX >> 2) + 1;    /* (2^(n/2-1))^2 = 2^(n-2) */
do	{
	if( x >= x2 )    /* x == x2 pode tomar qualquer ramo */
		{
		x2 = x2 + r + (h >> 1);
		r = (r+h) >> 1;
		}
	else
		{
		x2 = x2 - r + (h >> 1);
		r = (r-h) >> 1;
		}
	}
	while( (h>>=2) != 0 );
if( x <= (x2-r) )
	r--;
else if( x > (x2+r) )
	r++;

Comparação de Velocidade

Vamos analisar uma versão optimizada do código em assembly Intel APX/IA-32 (assumindo EAX=x e EBX=r, ambos inteiros sem sinal): 

	; …
	;                                                       ciclos i486:
	mov   EBX, 40000000h      ;  EBX = (2^(n/2-1)) << (n/2-1)          1
	mov   ECX, 20000000h      ;  ECX = (2^(n/2-1)) << (n/2-2)          1
	mov   EDX, 40000000h      ;  EDX = 2^(n-2)                         1
loop1:
	cmp   EAX, EDX            ;                                        1
	jb    bellow              ;                                      3/1
	;
	add   EDX, EBX            ;                                        1
	add   EBX, ECX            ;                                        1
	shr   EBX, 1              ;                                        2
	shr   ECX, 2              ;                                        2
	lea   EDX, [EDX+2*ECX]    ;                                        1
	jnz   loop1               ;                                      3/1
	jmp   done                ;                                      3/1
bellow:
	sub   EDX, EBX            ;                                        1
	sub   EBX, ECX            ;                                        1
	shr   EBX, 1              ;                                        2
	shr   ECX, 2              ;                                        2
	lea   EDX, [EDX+2*ECX]    ;                                        1
	jnz   loop1               ;                                      3/1
done:
	add   EDX, 2              ;  1 pelo ciclo, 1 para o "carry"        1
	sub   EDX, EBX            ;  EDX = limite inferior de EBX^2        1
	lea   ECX, [EDX+2*EBX]    ;  ECX = limite superior de EBX^2, +1    1
	cmp   EAX, EDX            ;                                        1
	sbb   EBX, 0              ;                                        1
	cmp   EAX, ECX            ;                                        1
	sbb   EBX, -1             ;                                        1

Isto é código de 32 bits, tal como se pode ver. Para converter este código para 16 bits, deveria:

  1. Substituir todos os registos de 32 bits por registos de 16 bits (i.e., remover o "E");
  2. Simplesmente remover os quatro zeros (0) mais à direita das primeiras três operações;
  3. Substituir os dois pares shr ECX, 2 / lea EDX, [EDX+2*ECX] por shr CX, 1 / add DX, CX / shr CX, 1 uma vez que uma deslocação por 2 e esse tipo de operação lea não são suportadas;
  4. Substituir o add EDX, 2 logo após a etiqueta (label) "feito" por inc DX;
  5. Substituir sub EDX, EBX / lea ECX, [EDX+2*EBX] logo após esse add por mov CX, DX / sub DX, BX / add CX, BX, uma vez que esse tipo de operação lea não é suportado.

O número de ciclos de relógio para este código assembly é de

Melhor caso:
3+12(n/2-1)-2+3+7 = 6n-1
Pior caso:
3+14(n/2-2)+12-2+3+7 = 7n-5

o que significa que este código pode calcular uma raiz quadrada de 32 bits (n=32) em 191~219 ciclos de relógio, e uma raiz quadrada de 16 bits (n=16) em 95~107 ciclos de relógio, o que é 14% a 30% mais rápido do que hard-sqrt!

Versão em Hardware

Se implementado em hardware, este algoritmo poderia ser paralelizado, o que significa ainda maior performance. Ambas as metades do ciclo seriam calculadas em paralelo, e então o resultado adequado seleccionado por um multiplexer a partir do resultado da comparação inicial. Deslocações não demoram tempo a executar já que podem ser implementadas ligando os bits de saída de um circuito, aos bits de entrada do registo, já deslocados.

Desde que haja tempo suficiente num ciclo de relógio para a propagação dos sinais eléctricos desde um circuito de adição / subtracção, um multiplexer e outro circuito de adição, o ciclo pode ser implementado em apenas 1 ciclo de relógio. Se não, usar multiplexers sincronizados com o relógio e fazendo algumas alterações, deveria permiti-lo implementar o ciclo em 2 ciclos de relógio, pelo menos.

Vejamos um diagrama do primeiro tipo de circuito:

[diagrama de circuito]

Este diagrama exibe apenas a parte do circuito da raiz quadrada que corre o ciclo principal — isto não inclui o código de prólogo ou epílogo. Entradas que têm uma nota de >>1 ou >>2 têm os bits do seu registo ligados de forma a implementar o deslocamento correspondente antes de entrar nesse circuito. Por exemplo se ligar o bit de saída u ao bit de entrada u-1 para todos os bits nessa palavra, e ligar o bit de entrada mais significativo a 0, está a implementar um deslocamento à direita por 1 entre os dois circuitos.

Todos os pinos de relógio de todos os registos (clk) devem ser ligado ao mesmo relógio comum. O ciclo termina (e EBX tem o resultado da raiz quase pronto), quando o pino de saída zero está a 1, indicando ECX=0. Os dois multiplexers devolvem o valor a se sel=0 (falso), ou b se sel=1 (verdadeiro). Finalmente, o circuito de adição que implementa EBX+ECX e o circuito de adição que se encontra mais à direita, podem ambos ser substituídos por circuitos OR de n bits.

Note que apesar do algoritmo ser altamente optimizado, o circuito que implementa o algoritmo não o é. Por exemplo, circuitos que executam quer uma adição, quer uma subtracção com base num sinal sel poderiam eliminar os dois pares ADD/SUB/MUX, substituindo cada por um único circuito ADD+SUB. De facto, o par ADD/SUB/MUX que implementa EBX±ECX pode ser substituído por um OR com ECX e um XOR especializado com o sinal sel. Um circuito de comparação especializado irá devolver o resultado EAX<EDX? muito mais rapidamente do que um circuito SUB, e assim por diante. O diagrama aqui exibido é uma ilustração de uma possibilidade, não uma optimização.

Virgula Flutuante

Isto estende-se facilmente a virgula flutuante (números reais). Números de virgula flutuante IEEE são guardados no formato

`+-"mantissa" * 2^"xpon",` `0 <= "mantissa" < 2,`
`|"xpon"| in NN_0`

onde xpon é o expoente e mantissa... a mantissa! Se o expoente tiver um valor especial chamado o valor "denormal", a mantissa deverá estar entre 0 (inclusive) e 1 (exclusive). Se o expoente tiver qualquer outro valor (um valor "normal"), a mantissa deverá estar entre 1 (inclusive) e 2 (exclusive). Isto significa que o bit mais significativo da mantissa pode ser determinado pelo expoente, pelo que a maior parte dos formatos de virgula flutuante não codificam este bit.

A mantissa pode ser transformada num inteiro multiplicando-a por uma potência de 2, o que é o mesmo de subtrair o expoente dessa potência de xpon. Isto reduz o problema ao cálculo da raiz quadrada de um inteiro (já sabemos fazer isto), e dividir xpon por 2 (trivial).

`sqrt("mantissa" * 2^"xpon") = sqrt("mantissa") * 2^("xpon"/2)`

De forma a devolver n bits significativos de dígitos, precisamos de trabalhar com 2n bits na mantissa, e ter o bit mais significativo da mesma a 1, pelo que precisamos de multiplicá-la por ainda outra potência de 2. Esta e outras pequenas coisas à volta desta questão (voltar a mantissa a um valor entre 0 e 2, lidar com um valor impar de xpon antes da sua divisão por 2, lidar com um número impar de bits na mantissa — alguns registos são inicializados a

`2^(n/2)`

e valores similares, que não são inteiros se n for impar! — etc.), são o que complica este processo.

Uma vez que está para além do objectivo deste documento, eu escrevi um programa de teste em ANSI C para descrever e testar dois algoritmos ligeiramente diferentes que calculam uma raiz quadrada em virgula flutuante (estes algoritmos devem ser idealmente implementados em hardware, e o programa de teste foi desenhado para facilitar a compreensão de como desenhar os mesmos nesse meio). Poderá fazer o download deste programa (código fonte em C, e versão em Inglês apenas) ao final desta secção.

O primeiro destes algoritmos é ligeiramente mais rápido e leva menos circuitos a implementar (o que significa que é mais barato), mas devolve valores que não são precisos, mas têm um pequeno (mas quantificável) erro. Ele deverá ser mais do que adequado para aplicações gráficas, por exemplo. Tal algoritmo pode ser implementado em hardware usando o seguinte diagrama:

[diagrama de circuito]

Os mesmos comentários que explicaram a implementação hardware do algoritmo de inteiros também se aplicam a este diagrama. Aqui, nomes de circuitos que terminam em *2 referem-se a circuitos com (ou que operam em) 2n bits em vez de apenas n bits. Notas de hi significam os n bits mais significativos de tais circuitos, enquanto low significam os n bits menos significativos. Por exemplo, uma das entradas ao circuito de adição mais à direita são os n bits mais significativos de ECX*2 deslocados para a direita por 1. O ciclo termina quando os n bits mais significativos de EBX*2 (sinal hi_zero) são 0. low(EBX*2) terá então a mantissa do resultado e xpon o seu expoente (que foi devidamente modificado antes do ciclo se iniciar). Para compreender a necessidade do contador para xpon, leia os comentários do código. A saída do circuito comparador SUB não necessita de ser complementada se trocar as entradas a e b dos multiplexers, mas o diagrama foi deixado como está de forma a que a única alteração estrutural ao diagrama do algoritmo inteiro seja nesta comparação.

O segundo algoritmo leva mais circuitos, mas não é mais complexo. Irá devolver valores exactos para as raízes quadradas.

[diagrama de circuito]

Tal como pode ver, agora quase todos os circuitos são de 2n bits, o que o torna quase duas vezes mais caro do que o algoritmo inteiro e acerca de 30% mais caro do que primeiro algoritmo de virgula flutuante já descrito. É no entanto o algoritmo ideal a implementar em FPUs (unidades de virgula flutuante, ou floating-point units, em Inglês), devido à sua precisão. O ciclo termina quando ECX*2 é 0 ou 1 (sinal zero_one), i.e., quando todos os seus bits à excepção do bit menos significativo são 0.

O programa de teste pode ser compilado para testar qualquer destes dois algoritmos em quase qualquer sistema compatível com ANSI C, e com quase qualquer tipo de formato de virgula flutuante (por favor leia a secção sobre restrições ao início do código). Poderá ser compilado em várias formas para testar uma variedade de formatos de virgula flutuante e variantes dos algoritmos. Pode testar valores pontuais (exibindo o progresso dos algoritmos) ou gamas inteiras de valores para uma prova empírica de precisão. Deve examinar o código da função gotest_real() para conhecer as implementações de ambos os algoritmos e seus comentários.

Por favor leia os comentários ao início do código para mais detalhes e conclusões (por exemplo, para números de virgula flutuante IEEE de 32 bits, o algoritmo de hardware seria aproximadamente 30% mais rápido do que um Intel Itanium de ultima geração (em Agosto de 2000)).

  Versão 1.00
ftest.zip (21kb)

Raiz de Qualquer Índice

Quando se usa o algoritmo de pesquisa binária para uma raiz de índice i, os valores das raízes que se tomam a cada iteração (no primeiro exemplo, 128 seguido de 64 ou 192, seguido de 32 ou ..., e assim por diante) não são idênticos. O primeiro valor para o qual se irá testar a raiz (iremos chamar-lhe m) é metade da raiz máxima que podemos gerar usando números de n bits, ou seja

`m = root(i) (2^n-1) / 2`

o que é muito próximo de

`m = 1/2 * root(i) (2^n)`

`m = 2^(n/i-1)`

Se este valor não for inteiro, necessitará de usar um m double em vez de int. r e x2 no código C descrito abaixo seguem o mesmo tipo de m. Para n=16 e i=2 (o exemplo que temos estado a usar), m=128, tal como esperado.

Este valor ainda é dividido por 2 a cada iteração do ciclo, mas ao contrário do que acontecia com i=2, onde m era sempre divisível por 2, já não podemos garantir isso pelo que precisamos de usar m/=2 e m double em vez de m>>=1 nesses casos.

Actualizar o valor de x2 dentro do ciclo C também já não é idêntico, tal como poderia esperar. O problema não é complexo, no entanto. Se chamar a XA(r,m)

`(r+m)^i - r^i`

e a XS(r,m)

`(r-m)^i - r^i`

então pode adicionar estas fórmulas a x2 ao adicionar ou subtrair m de r, respectivamente. Por exemplo, para i=2, XA(r,m)=m*(2*r+m), e para i=3, XA(r,m)=m*(3*r*(r+m)+m*m).

Finalmente, precisamos de uma condição para terminar o ciclo. Poderíamos usar m<0.001 ou outro valor baixo, mas precisaríamos ajustar este valor para cada índice. Também poderíamos usar XA(r,m)<0.5 mas essa condição iria forçar a demasiadas iterações. A solução encontrada foi terminar o ciclo quando a) o "ramo" do ciclo (adição/subtracção) mudou, e b) r arredondado ao inteiro mais próximo não mudou com o ramo. Isto significa que a raiz está entre os dois quadrados previamente testados, e é igual a r arredondado ao inteiro mais próximo. Isto permite-nos eliminar o código de epílogo (o que não é favorável em termos de velocidade, mas simplifica o código), mas tem uma desvantagem — os valores mínimo e máximo para números de n bits vão causar um ciclo infinito, já que o ciclo vai continuamente aproximar-se deles, mas nunca ultrapassá-los. Eles deveriam ser tratados como casos especiais.

Com isto, podemos agora criar o seguinte código em C: 

#include <math.h>
#include <limits.h>
/* … */
#define  n  ((int) floor( log(ULONG_MAX) / log(2) ))
#define  i  /* defina i como sendo o indice da raiz */
/* … */
register unsigned long int x, prev_r;
register int prev_side;
register double r, h, x2;

/* … */
r  = pow( 2, (double) n/i-1 );    /* 2^(n/i-1) */
x2 = (ULONG_MAX >> i) + 1;        /* r^i = 2^(n-i) */
/* nota: recodifique a linha "r=" acima de forma a que os calculos
   sejam feitos ou por voce^ ou pelo compilador durante a compilacao.
   compiladores inteligentes *podem* ja' fazer isso */
h = r;
prev_side = 0;    /* nenhum; logo, prev_r nao precisa inicializacao */
for(;;)    /* ciclo infinito */
	{
	h /= 2;
	if( x >= x2 )    /* x == x2 pode tomar qualquer ramo */
		{
		if( prev_side == -1 )
			{
			if( prev_r == ((int) r+.5) )
				break;
			prev_side = +1;
			}
		prev_r = ((int) r+.5);
		x2 += XA(r,h);
		r += h;
		}
	else
		{
		if( prev_side == +1 )
			{
			if( prev_r == ((int) r+.5) )
				break;
			prev_side = -1;
			}
		prev_r = ((int) r+.5);
		x2 += XS(r,h);
		r -= h;
		}
	}
r = (int) r+.5;

Este código irá funcionar apenas na gama 0<x<ULONG_MAX (i.e., números positivos, com excepção dos casos especiais 0 e ULONG_MAX). Qual a eficiência deste código? Depende apenas da sua capacidade de optimizar código. Por favor veja a secção de Raiz de Qualquer Índice do algoritmo sequencial na Parte I para mais detalhes e exemplos. Não irei aprofundar esta questão.

Trabalho a Desenvolver

A raiz quadrada é frequentemente usada em fórmulas (especialmente fórmulas 3D) multiplicada por uma constante ou outro valor. Uma vez que o algoritmo da pesquisa binária itera por cada bit do resultado e logo é semelhante à funcionalidade de um circuito de multiplicação simples, poderia-se alterar o primeiro (ou segundo) diagrama de virgula flutuante aqui exibido para que uma multiplicação fosse executada em paralelo com o cálculo da raiz quadrada, e fosse devolvido um resultado inteiro (ou de virgula flutuante, se modelar o circuito para tal) com muito pouca perda de precisão (se alguma).

Isto é relativamente simples de executar se pensar nisso: enquanto EBX está a ser construído em qualquer variante do algoritmo, o bit representado por ECX irá sempre tornar-se 1 em EBX, independentemente se lhe adicionou ou subtraiu ECX, mas o bit anterior a esse (e apenas esse bit) poderá ser modificado (se foi usada uma subtracção) de 1 para 0. A multiplicação paralela poderá então seguir esse bit anterior em cada ciclo, em vez do bit em ECX. Não irei aprofundar este assunto, no entanto — o processo é linear e pode dar origem a circuitos gráficos 3D muito simplificados e acelerados.

Outra direcção de trabalho é no desenvolvimento uma "Pesquisa Estimada". Esta pesquisa foi conceptualizada devido aos primeiros dois parágrafos que descrevem a pesquisa binária.

Quando você pesquisa um dicionário, você não o abre logo ao meio, mas na realidade faz uma estimativa de onde pensa que a palavra poderá estar. E depois, pela comparação que faz com a palavra que lá encontrou, você não corta cegamente uma das duas secções ao meio, mas na realidade move-se para dentro delas uma certa distância — mais uma vez uma estimativa.

Mas com que precisão é que consegue estimar raízes? Bem, como já vimos, "comer" metade dos dígitos menos significativos de um número (ou a parte inteira dessa metade arredondada para baixo se o número de dígitos significativos é impar) devolve uma aproximação bastante boa. Isto vem de

`sqrt(2^n) = 2^(n/2)`

Assim, se encontrar um n tal que

`2^(n-1) <= x < 2^n`

`n = "int"(log_2 x)+1`

onde o logaritmo é arredondado para baixo, então

`x / 2^(n/2)`

seria um aproximação bastante boa (e rápida). O erro é dado por

`sqrt(x) - x / 2^(n/2)`

que tem o seu valor mínimo a quase

`sqrt(2^n) - 2^n / 2^(n/2) = 0`

tal como esperado, e o seu máximo a

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

o que para palavras de 16 bits que têm o n máximo a 16 (isto é um significado ligeiramente diferente de n do que aquele que lhe temos dado até agora), o erro máximo é de 54 (>21%); para palavras de 32 bits, o n máximo é 32 e o erro máximo 13573 (<21%).

A seguinte amostra de código assembly Intel APX/IA-32 (x86) devolve a aproximação em EBX, deixando EAX sem modificar, mas destruindo ECX. 

	; …
	;                                       i486 cycles:
	bsr   ECX, EAX    ;  ECX = n-1                 6-103
	inc   ECX         ;  ECX = n                       1
	shr   ECX, 1      ;                                2
	mov   EBX, EAX    ;  2^(n-1) <= EAX,EBX < 2^n      1
	shr   EBX, CL     ;  EBX = EBX / 2^(n>>1)          2
	; …

Infelizmente, você não pode simplesmente adicionar este código ao topo de um algoritmo sequencial tal como está. A aproximação impediu o algoritmo sequencial de reduzir EAX tal como normalmente o faria. Adicioná-lo ao algoritmo de pesquisa binária também não é boa ideia porque o erro máximo que tem significa que apenas uma iteração do algoritmo seria evitada, e demasiado tempo é gasto no cálculo da aproximação. Por outro lado, a única forma de encontrar o quadrado desta aproximação (de forma a poder ser usado em qualquer destes algoritmos) é fazendo uma multiplicação o que não é tão rápido quanto isso. Mais trabalho nesta questão é assim deixado ao leitor.

Conclusão

Em termos de velocidade, a pesquisa binária é um dos melhores (ou o melhor) algoritmo a usar por software em CPUs que apenas tenham operações inteiras, ou para se facilmente implementar em hardware. Em termos de espaço, a pesquisa sequencial crescente descrita na Parte I não tem igual.

Talvez estes não sejam os melhores algoritmos existentes, e eu não serei (possivelmente) o primeiro a desenvolvê-los. Desconheço o seu uso onde quer que seja, no entanto, ou documentos ou livros que descrevam algoritmos tais como os que eu aqui descrevo, mas não efectuei pesquisa extensiva sobre esta questão.

Quanto à minha contribuição: espero que goste.

Fotografia pela NASA no Unsplash