Social Icons

^^

segunda-feira, 15 de agosto de 2011

Delphi com assembly (Parte 6)

 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;

Nenhum comentário:

Postar um comentário

Popular Posts

 

Seguidores

Hora exata:

Total de visualizações de página