Como funciona a aritmética fundamental em ponto flutuante

Quando falamos em aritmética fundamental usando valores inteiros (adição, subtração, multiplicação e divisão) acredito que não há grandes dúvidas, quando os argumentos estão expressos em formato binário. Mas, e quanto à codificação em ponto flutuante que citei em artigos anteriores?

Adição e subtração

O problema que enfrentamos é que, além dos bits significativos, temos o escalonamento (ou o “deslocamento” do ponto). Na adição precisamos que ambos os valores tenham a mesma escala para poder adicioná-los. Um exemplo, em decimal, é a tentativa de adicionar 1.0\cdot10^0 ao valor 1.0\cdot10^2. Não podemos, simplesmente, somar 1.0 com 1.0! Temos que tornar as escalas (10^0 e 10^2) iguais e, para isso, adaptamos o menor valor:

\displaystyle 1.0\cdot10^2 + 0.01\cdot10^2=1.01\cdot10^0

Aqui, 1.0\cdot10^0 foi transformado em 0.01\cdot10^2. O motivo de adequarmos o menor valor é que os bits do menor valor serão deslocados para a direita e a adição tenderá a obter uma soma normalizada. Troque agora o valor normalizado por um valor inteiro binário (com o bit 23, implícito, setado) e a escala na base 2, ao invés de 10… Lembre-se que, no artigo anterior, demonstrei que o valor contido na “fração” da estrutura de um float é o valor inteiro do numerador, cujo denominador é sempre o mesmo, 2^{23}.

Embora esteja mostrando o funcionamento do ponto de vista de um valor fracionário, o coprocessador realiza a adição com o valor inteiro. Usando o programinha shwflt, do artigo anterior, temos:

$ ./shwflt 1e0
[Normal] 1 -> s=+, E=0, f=0 (0b1.00000000000000000000000)
$ ./shwflt 1e2
[Normal] 100 -> s=+, E=6, f=0x480000 (0b1.10010000000000000000000)

Consideremos a=1.0\cdot10^2 e b=1.0\cdot10^0. Com o procedimento descrito acima, temos que a>b e, portanto, b deve ser adequado. Em binário, considerando apenas as frações e o bit implícito, temos, em binário, a=0b110010000000000000000000 e b=0b100000000000000000000000, mas b deve ser deslocado 6 bits para a direita para que os expoentes sejam equalizados (E_a=E_b) e, portanto, b=0b000000100000000000000000. A adição das frações, então, dá a+b=0b110010100000000000000000 e o resultado final, na estrutura do float será f=0b10010100000000000000000,\,E=6, como pode ser visto aqui:

$ ./shwflt 1.01e2
[Normal] 101 -> s=+, E=6, f=0x4a0000 (0b1.10010100000000000000000)

Existe um detalhe interessante: Se a diferença entre os expoentes E_a e E_b for maior que 24, significa que o menor valor pode ser seguramente desconsiderado na adição. O motivo, óbvio, é que ele deverá ser deslocado 24 bits para a direita, fazendo com que todos os seus bits significativos “desapareçam”. Assim, a adição só é feita se E_a - E_b < 24 onde E_a > E_b.

A subtração segue as mesmas regras, mas a ordem com a qual os operandos aparecem é importante.

Renormalização

Consideremos o caso de adicionarmos o valor fracionário, binário, 0b1.0 a ele mesmo. Teremos que obter o valor 2.0 ou 0b10.0. No entanto o bit MSB é implícito e precisa estar isolado da parte fracionária. Por isso, na adição, precisamos ter espaço suficiente para o armazenamento para um bit extra à esquerda, no valor inteiro. Ou seja, o resultado por ter 25 bits de tamanho: Assim, se R=0x800000+0x800000=0x1000000, o valor de R precisa ser deslocado 1 bit para a direita (R\text{ shr }1 e o expoente E é incrementado. Vejamos:

$ ./shwflt 2
[Normal] 2 -> s=+, E=1, f=0 (0b1.00000000000000000000000)

Repare que f=0, do mesmo jeito que acontece com o valor 1.0, mas o expoente E=1 por causa da renormalização.

Da mesma forma, isso também pode ocorrer com a subtração. O MSB sempre tem que ser igual a 0b1, a não ser que obtenhamos um valor subnormal. Se o MSB obtido for zero o valor final será deslocado para a esquerda tantos bits quantos forem necessários para que ele seja igual a 0b1, mas o expoente E, do resultado, será subtraído dessa mesma quantidade de bits, normalizando o valor. Se o coprocessador perceber que E será menor que -126, então o valor é subnormal e a renormalização não é feita.

Multiplicação e divisão

Diferente da adição (e da subtração), a multiplicação não precisa sofrer adequações das frações. Isso fica claro quando você percebe que basta somar os dois expoentes de mesma base para obter a escala correta e, as frações, só precisam ser multiplicadas:

\displaystyle (f_a\cdot2^{E_a})\cdot(f_b\cdot2^{E_b}) = (f_a\cdot f_b)\cdot2^{E_a + E_b}

O detalhe é que, ao obtermos a fração com o bit implícito, teremos um valor de 24 bits de tamanho que, ao multiplicar por outro valor de 24 bits de tamanho nos dará um valor de 48 bits de tamanho. Considere o caso de multiplicarmos 1.0 por ele mesmo. Esse valor é codificado, em hexadecimal, como 0x800000. Depois da multiplicação obtemos 0x400000000000. Não é óbvio, mas o resultado encontra-se 23 bits deslocado para a esquerda, ou seja, justamente os 23 bits da parte fracionária dos argumentos originais! A resposta é justamente \left[(a\cdot b)\text{ shr }23\right]\cdot2^{E_a+E_b}, onde a e b são os valores inteiros de 24 bits com o bit implícito!

A divisão, é claro, segue o sentido inverso. os bits significativos devem ser divididos e os expoentes subtraídos, mas diferente da multiplicação, o dividendo deve ser deslocado 23 bits para a esquerda antes que a divisão possa ser feita. Novamente, vamos tentar dividir 1.0 por ele mesmo… Se dividirmos 0x800000 por ele mesmo obteremos o quociente 1. Mas 1 é apenas o bit 0 da fração, ou seja, 2^{-23}. Mas, se deslocarmos a para esquerda em 23 bits, teremos \frac{0x400000000000}{0x800000}=0x800000. Assim, a divisão binária fica \frac{a\text{ shl }23}{b}\cdot2^{E_a-E_b}. Claro que temos que considerar o caso espacial onde $label b=0$.

Atenção que a renormalização também pode ser necessária depois de multiplicações e divisões. Tomemos o caso onde todos os bits significativos estejam setados e E=0. Se quisermos elevar esse valor ao quadrado, neste caso, a será 0xffffff e obteremos 0xfffffe000001. Ao deslocarmos 23 bits para a direita teremos 0x1fffffc, que tem 25 bits de tamanho! Isso deve ser deslocado em 1 bit para a direita e o valor de E incrementado. De fato, o programinha abaixo atesta esse fato:

#include <stdio.h>

union flt_u {
  float v;
  struct {
    unsigned int f:23;
    unsigned int e:8;
    unsigned int s:1;
  };
};

void main(void)
{
  /* Todos os bits significativos setados.
     E=0. 
     s=+ */
  union flt_u flt = { f:~0, e:127 };

  float b;

  b = flt.v * flt.v;

  printf("a = %.18f, a*a = %.18f, "
         "f[a*a]=%#x, E[a*a]=%d\n",
    flt.v,
    b,
    (*(union flt_u *)&b).f,
    (*(union flt_u *)&b).e-127);
}

Ao executar o código acima, temos:

a = 1.999999880790710449, a*a = 3.999999523162841797, f[a*a]=0x7ffffe, E[a*a]=1

Como os dois valores originais tẽm expoente zerado, o expoente final deveria ser zerado também (E_a+E_b=0+0=0, mas a renormalização nos deu um expoente unitário, sendo fácil perceber que a parte fracionária é resultante do shift para a direita do valor que expus anteriormente. Na divisão algo similar pode acontecer, mas no sentido contrário.

Underflows. overflows e NaNs

O que acontece se tivermos f_a com todos os bits setados e E_a=127 (o maior expoente possível em base 2), pegarmos esse valor máximo e adicionarmos 2^{104}? O valor 2^{104} é precisamente o bit 0 da estrutura do float com um expoente 127. O que ocorrerá que o coprocessador fará o cálculo, mas obterá um valor “estranho” com expoente 128 que, por definição do padrão, corresponde a ∞ ou a NaN. No caso da adição, ao ocorrer um overflow termos o resultado ±∞, dependendo do sinal. Mas, divisões podem causar resultados NaN. É o caso da divisão \frac{0}{0}.

Ao obter, no resultado, um expoente maior que 127 ou menor que -127 (para valores subnormais) o coprocessador decide o que vai deixar na fração… Se o valor for “infinito”, a fração será zerada. Se for um NaN, geralmente temos algum “lixo” armazenado por lá… O bit 22 será 0 se tivermos um qNaN e 1 se tivermos um sNaN, mas não podemos usar os demais bits para qualquer coisa, senão forçarmos um NaN (para o caso do qNaN).

Parece “ponto fixo”, não é?

Não sei se escrevi algo sobre ponto fixo neste blog (espero que sim, me avisem caso contrário, ok?), mas esse esquema de usarmos os bits significativos como se fossem valores inteiros lembra muito a maneira como ponto fixo funciona, com uma grande vantagem: Se dividirmos os 32 bits em partes iguais entre o componente inteiro e o fracionário poderemos, em ponto fixo, expressar o valor mínimo de 2^{-16} e o máximo de \frac{2^{32} - 1}{2^{16}}. Ponto flutuante, do mesmo tamanho, em contrapartida, nos garante a faixa entre o valor mínimo subnormal de 2^{-149}\approx1.4\cdot10^{-45} até o máximo de \frac{2^{25} - 1}{2^{23}}\cdot2^{127}\approx6.8\cdot10^{38}. Uma faixa bem maior e com maior precisão usando menos bits significativos (apenas 24).

Note que, para ser exato, os 32 bits de um ponto fixo deveriam ser divididos em 31 bits significativos e um de sinal. E, no neste caso, é comum usarmos o bit de sinal para conter valores inteiros em complemento 2, coisa que não é feita com ponto flutuante. Os bits significativos, neste último caso, são sempre armazenados de forma positiva. Fica a cargo dos algoritmos determinar se, por exemplo, estamos realizando operações com valores com sinais iguais ou diferentes, o que adiciona alguma complexidade (pouca!). E, é claro, o padrão IEEE 754 ainda fornece alguns valores especiais, o que não é diretamente possível com ponto fixo.

Anúncios