V¡®u§ 849 Posted January 7, 2012 Aritmética Inteira de 128 bits (3) Nessa terceira e última parte da série sobre inteiros de 128-bit (que chamamos de inteiros gigantes, inteiros enormes ou inteiros largos) vamos finalmente lidar com a aritimética propriamente dita- as quatro operações fundamenteais (adição, subtração, multiplicação e divisão). Antes de iniciar, gostaria de avisar que as rotinas introduzidas nas duas partes anteriores foram corrigidas e também otimizadas. Ainda não pude testá-las tanto quanto gostarias. Se você encontrar erros ou tiver algum comentário sobre os fontes, por favor me envie um e-mail. Adição ====== Como somamos dois números, cada um formado de quatro inteiros de 32-bit? Bem, é bem fácil. Simplesmente somamos como faríamos com dois números de quatro dígitos decimais (como 3597 e 0015, por exemplo), exceto que cada dígito aqui pode ter aproximadamente 4 bilhões (2^32) de valores distintos ao invés de apenas dez. O algoritmo é algo como: function AddWithCarry(x: Longint; y: Longint; var Carry: Boolean): Longint; forward; function HugeAdd(x: Hugeint; y: Hugeint): Hugeint; // Result := x + y; var Carry: Boolean; begin Carry := False; Result[0] := AddWithCarry(x[0], y[0], Carry); Result[1] := AddWithCarry(x[1], y[1], Carry); Result[2] := AddWithCarry(x[2], y[2], Carry); Result[3] := AddWithCarry(x[3], y[3], Carry); end; AddWithCarry é uma função fictícia que retorna um inteiro com os 32 bits de mais baixa ordem do resultado da adição dos dois argumentos, mais um se Carry (o terceiro argumento) for True. A rotina também altera o valor de Carry com True o False dependendo da adição ter gerado um excesso ou não (ou se o carry é 1 ou 0, se você prefere ver dessa forma). Na prática, essa função não precisa ser fictícia: function AddWithCarry(x: Longint; y: Longint; var Carry: Boolean): integer; asm // if Carry then CF := 1 else CF := 0; test byte ptr [ecx], -1 // efeito colateral: CF := 0; jz @@NoCarry stc // CF := 1; @@NoCarry: // Result := x + y + CF; CF := GeneratedCarry; adc eax, edx // Carry := CF; setc byte ptr [ecx] end; é mais eficiente codificar HugeAdd inteiramente em assembler: function HugeAdd(x: Hugeint; y: Hugeint): Hugeint; // Result := x + y; // Parameters: EAX = @x; EDX = @y; ECX = @Result asm push esi mov esi, [eax+_0_] // ESI := x[0]; add esi, [edx+_0_] // ESI := ESI + y[0]; mov [ecx+_0_], esi // Result[0] := ESI; mov esi, [eax+_1_] // ESI := x[1]; adc esi, [edx+_1_] // ESI := ESI + y[1] + Carry; mov [ecx+_1_], esi // Result[1] := ESI; mov esi, [eax+_2_] // ESI := x[2]; adc esi, [edx+_2_] // ESI := ESI + y[2] + Carry; mov [ecx+_2_], esi // Result[2] := ESI; mov esi, [eax+_3_] // ESI := x[3]; adc esi, [edx+_3_] // ESI := ESI + y[3] + Carry; mov [ecx+_3_], esi // Result[3] := ESI; pop esi end; Subtração ========= A subtração acontece de forma semelhante à adição mas, ao invés de gerar um excesso, a operação gera um "empréstimo" (também representado pelo flag Carry) se o minuendo (primeiro operando) é menor que o subtraendo (segundo operando): function SubtractWithBorrow(x: Longint; y: Longint; var Borrow: Boolean): Longint; forward; function HugeSub(x: Hugeint; y: Hugeint): Hugeint; // Result := x - y; var Borrow: Boolean; begin Borrow := False; Result[0] := SubtractWithBorrow(x[0], y[0], Borrow); Result[1] := SubtractWithBorrow(x[1], y[1], Borrow); Result[2] := SubtractWithBorrow(x[2], y[2], Borrow); Result[3] := SubtractWithBorrow(x[3], y[3], Borrow); end; function SubtractWithBorrow(x: Longint; y: Longint; var Borrow: Boolean): Longint; asm // if Borrow then CF := 1 else CF := 0; test byte ptr [ecx], -1 // Side-effect: CF := 0; jz @@NoBorrow stc // CF := 1; @@NoBorrow: // Result := x - y - CF; CF := NeededBorrow; sbb eax, edx // Borrow := CF; setc byte ptr [ecx] end; Você deve estar pronto para escrever uma versão de HugeSub inteiramente em assembler já que é análoga à HugeAdd; tudo que você deverá fazer é substituir ADD e ADC com SUB e SBB respectivamente. Número Oposto ============= Dado um número, essas implementações de HugeNeg retornam seu oposto (complemento de dois): function HugeNeg(x: Hugeint): Hugeint; begin // Result := (Not x) + 1; Result := HugeAdd(HugeNot(x), IntToHuge(1)); end; function HugeNeg(x: Hugeint): Hugeint; begin // Result := 0 - x; Result := HugeSub(IntToHuge(0), x); end; A segunda é mais simples e rápida pois involve uma operação apenas e, já que sabemos subtrair, podemos implementá-la em assembler: function HugeNeg(x: Hugeint): Hugeint; // Result := -x; // Parameters: EAX = @x; EDX = @Result asm // Result := 0 - x; push esi xor esi, esi mov ecx, [eax+_0_] // x[0] sub esi, ecx // 0 - x[0] mov ecx, 0 mov [edx+_0_], esi // Result[0] mov esi, [eax+_1_] // x[1] sbb ecx, esi // 0 - x[1] - Borrow mov esi, 0 mov [edx+_1_], ecx // Result[1] mov ecx, [eax+_2_] // x[2] sbb esi, ecx // 0 - x[2] - Borrow mov ecx, 0 mov [edx+_2_], esi // Result[2] mov esi, [eax+_3_] // x[3] sbb ecx, esi // 0 - x[3] - Borrow mov [edx+_3_], ecx // Result[3] pop esi end; Multiplicação ============= Uma forma de multiplicar números é através de um laço de adições: function HugeMul(x: Hugeint; y: Hugeint): Hugeint; begin SetZero(Result); while not HugeIsZero(y) do begin Result := HugeAdd(Result, x); HugeSub(y, 1) end; end; Em termos computacionais, esse algoritmo é bastante pobre. Por exemplo, se o valor de "y" for 4 milhões, o laço se repetiria 4 milhões de vezes! De qualquer forma, a idéia ainda é boa se pudéssemos de alguma forma acelerar o processo. Vamos então brincar um pouco com álgebra de bits: x * y = x * (y[3]*2^96 + y[2]*2^64 + y[1]*2^32 + y[0]*2^0) = (x*y[3])*2^96 + (x*y[2])*2^64 + (x*y[1])*2^32 + (x*y[0])*2^0 Agora nós reduzimos o problema de multiplicar dois inteiros gigantes ao de multiplicar um inteiro gigante por um inteiro de 32-bit. Quando multiplicamos o primeiro operando pelos quatro inteiros que compõem o segundo operando, então deslocamos os resultados parciais por 0, 32, 64 e 96 bits (multiplicação por 2^0, 2^32, 2^64 e 2^96) e finalmente somamos esses valores para obter o resultado final. function HugeMulInt(x: Hugeint; y: Longint): Hugeint; forward; function HugeMul(x: Hugeint; y: Hugeint): Hugeint; begin Result := HugeShl(HugeMulInt(x, y[3]), 96) + HugeShl(HugeMulInt(x, y[2]), 64) + HugeShl(HugeMulInt(x, y[1]), 32) + HugeMulInt(x, y[0]); end; Essa é a exata forma que multiplicamos números decimais quando usamos papel, exceto que aqui a base é 2^32 ao invés de dez. Vejamos como multiplicar um inteiro gigante por um inteiro: function MultiplyWithCarry(x: Longint; y: Longint; var Carry: Longint): Longint; forward; function HugeMulInt(x: Hugeint; y: Longint): Hugeint; // Result := x * y; var Carry: Longint; begin Carry := 0; Result[0] := MultiplyWithCarry(x[0], y, Carry); Result[1] := MultiplyWithCarry(x[1], y, Carry); Result[2] := MultiplyWithCarry(x[2], y, Carry); Result[3] := MultiplyWithCarry(x[3], y, Carry); end; function MultiplyWithCarry(x: Longint; y: Longint; var Carry: Longint): integer; // Result := LoDWord(x * y + Carry); // Carry := HiDWord(x * y + Carry); asm // EDX:EAX := EAX * EDX; // x * y mul edx // Inc(EDX:EAX, Carry); add eax, [ecx] adc edx, 0 // Carry := EDX; // High order 32 bits of the result mov [ecx], edx; end; MultiplyWithCarry é muito semelhante a AddWithCarry mas realiza uma multiplicação ao invés de uma adição, e gera um excesso de 32-bit ao invés de um bit apenas (a multiplicação de dois valores de 32-bit gera um valor de 64-bit enquanto sua soma gera um valor de 33-bit). MultiplyWithCarry primeiramente realiza a multiplicação sem sinal de "x" (EAX) por "y" (EDX), usando o opcode MUL. O resultado de 64-bit é um inteiro sem sinal contido em EDX:EAX, ao qual a função soma o valor do parâmetro Carry. A função retorna os 32 bits mais baixos do resultado final (localizado em EAX) e os 32 bits mais altos (EDX) constituem o excesso para a próxima multiplicação, que é armazenado no parâmetro Carry (passado por referência). Uma implementação em assembler de HugeMul e HugeMulInt encontra-se no código em anexo. Por razões de simplificação, nos exemplos acima, as funções consideram que os números não possuem sinais, mas o código anexo considera os sinais. O código de HugeMul também não chama HugeMulInt ou HugeShl e está bastante otimizado. Ao invés de considerar um inteiro gigante como quatro inteiros de 32-bit, podemos considerá-lo como 128 inteiros de 1-bit multiplicados por 128 potências de 2: bit127 * 2^127 + bit126 * 2^126 + ... + bit1 * 2^1 + bit0 * 2^0 Como cada bit só pode ser 0 ou 1, o algoritmo acima pode ser muito simplificado: function HugeMul(x: Hugeint; y: Hugeint): Hugeint; // Result := x * y; var i: Longint; begin SetZero(Result); for i := 0 to 127 do if BitTest(y, i) then Result := HugeAdd(Result, HugeShl(x, i)); end; A idéia é somar diferentes potências de 2 de "x", dependendo essas potências nos bits de "y". Por exemplo, se "y" for 20, os bits 5 e 3 estariam ligados (20 em decimal é 10100 em binário), de modo que apenas duas adições são realizadas e o resultado é HugeShl(x, 3) mais HugeShl(x, 5). Esse algoritmo pode ser codificado de forma bastante eficiente em assembler mas o primeiro algoritmo será ainda mais rápido. A razão de mostrá-lo é que facilitará o entendimento do algoritmo que usaremos para divisões. Divisão ======= Vamos primeiro ver o caso da divisão de um Hugeint por um inteiro de 32-bit, que deve ser fácil de entender: function DivideWithRemainder(x: Longint; y: Longint; var Remainder: Longint): Longint; forward; function HugeDivInt(x: Hugeint; y: Longint): Hugeint; // Result := x div y; var Remainder: Longint; begin Remainder := 0; Result[0] := DivideWithRemainder(x[3], y, Remainder); Result[1] := DivideWithRemainder(x[2], y, Remainder); Result[2] := DivideWithRemainder(x[1], y, Remainder); Result[3] := DivideWithRemainder(x[0], y, Remainder); asm mov edx, Remainder end; end; function DivideWithRemainder(x: Longint; y: Longint; var Remainder: Longint): Longint; // Result := Remainder:x div y; // Remainder := Remainder:x mod y; asm push esi mov esi, edx // y mov edx, [ecx] // Remainder // EAX := EDX:EAX div ESI; // EDX := EDX:EAX mod ESI; div esi // Remainder := EDX; mov [ecx], edx; pop esi end; HugeDivInt deixa o resto da divisão em EDX, de modo que possa ser usado em uma função que retorne o resto da divisão: function HugeModInt(dividend: Hugeint; divisor: Longint): Longint; // Result := dividend mod divisor; // Parameters: EAX = @dividend; EDX = @divisor; asm sub esp, TYPE(Hugeint) // Abrir espaço na pilha para um Hugeint mov ecx, esp // com o resultado da divisão call HugeDivInt // Realizar a divisão add esp, TYPE(Hugeint) // Restaurar o ponteiro da pilha mov eax, edx // Result := Remainder; // o resto foi deixado em EDX end; Para o caso de dois inteiros gigantes, temos que pensar no algoritmo que usaríamos para dividir um número de quatro dígitos decimais no papel. Mas esse procedimento é complexo e lento demais pois envolve além das divisões, multiplicações, subtrações, antecipação de etapas, etc. Há alternativa possível para o algoritmo? Sim: function HugeDiv(dividend: Hugeint; divisor: Hugeint): Hugeint; // Result := dividend div divisor; begin if HugeIsZero(divisor) then raise EDivByZero.CreateRes(@sDivByZero); Result := 0; while HugeCmp(dividend, divisor) >= 0 do begin dividend := HugeSub(dividend, divisor); Result := HugeAdd(Result, IntToHuge(1)); end; end; Claro, o algoritmo acaba sendo extremamente lento (se dividirmos 12 milhões por 3, o laço executará 4 milhões de vezes) mas podemos acelerá-lo se subtrarimos do dividendo o divisor multiplicado por diferentes potências de 2, da maior para a menor, definindo o bit correspondente do resultado toda vez que efetuarmos a subtração (o bit na posição da potência de 2 que foi usado). é o inverso do que foi feito no caso da multiplicação acima. O processo de divisão seria reduzido ao máximo de 128 subtrações. No exemplo a seguir, o dividendo é 20 (10100 em binário) e o divisor é 3 (11 em binário): 10100 - 11 * 2^2 = 10100 - 1100 = 1000 Result := 100 1000 - 11 * 2^1 = 1000 - 110 = 10 Result := 110 Inicialmente, 11 * 2^2 é o maior valor que é menor ou igual ao dividendo então subtraímos esse valor do dividendo e definimos o bit 2 do resultado pois subtraímos o divisor multiplicado por 2^2. Até o momento, o resto da divisão é 8 (1000 em binário) e 11 * 2^1 é o maior valor que é menor ou igual ao resto. Então subtraímos esse valor do resto e definimos o bit 1 do resultado pois subtraímos o divisor multiplicado por 2^1. O resto agora é 2 (10 em binátio) e como o divisor é maior que esse valor, a divisão é interrompida. O resto da operação seria 2 (10 em binário) e, como os bits 2 e 1 do resultado estão marcados, o resultado é 110 em binário, ou 6 em notação decimal. function HugeDiv(dividend: Hugeint; divisor: Hugeint): Hugeint; var _r_: Hugeint; // remainder _d_: Hugeint; // divisor _q_: Hugeint; // quotient BitPosR, BitPosD, count: integer; begin _r_ := dividend; _d_ := divisor; HugeSetZero(_q_); BitPosD := HugeBitScanReverse(_d_); if BitPosD = -1 then RaiseDivByZero; BitPosR := HugeBitScanReverse(_r_); count := BitPosD - BitPosR; if count > 0 then _d_ := HugeShl(_d_, count); repeat if HugeCmp(_d_, _r_) <= 0 then begin _r_ := HugeSub(_r_, _d_); HugeBitSet(_q_, count); end; _d_ := HugeShr(_d_, 1); dec(count); until count < 0; Result := _q_; asm lea edx, _r_ end; end; HugeBitScanReverse é a função que retorna a posição do primeiro bit não zero, realizando a busca a partir do bit 127 até o bit 0. Se todos os bits são zero, o resultado é -1. Usamos HugeBitScanReverse para determinar a primeira potência de dois que devemos multiplicar o divisor para iniciar a divisão. A implementação em assembler de HugeDiv no código anexo suporta inteiros com sinais. é apenas uma primeira aproximação que pode ainda ser muito melhorada. A função deixa em EDX o endereço do resto de modo que ele possa ser usado por uma função que retorna o módulo da divisão: function HugeMod(dividend: Hugeint; divisor: Hugeint): Hugeint; // Result := dividend Mod divisor; // Parameters: EAX = @dividend; EDX = @divisor; ECX = @Result asm push ecx // @Result call HugeDiv // EDX := @remainder; pop eax // EAX := @Result; call HugeMov // EAX^ := EDX^; end; Share this post Link to post Share on other sites
