Valores exatos em ponto flutuante

Eu já escrevi bastante, por aqui, a respeito de ponto flutuante. Fiz isso porque um monte de gente acha que usar float ou double é a solução ideal para garantir a “precisão” de cálculos com “valores quebrados”, mesmo não entendendo o que a palavra “precisão” significa nesse contexto… Se você consultar o blog verá, pelo menos, dez artigos à respeito. Para aumentar um cadinho a lista vou citar uma vantagem, não um problema.

Topei, recentemente, com um projeto onde o desenvolvedor substitui a estrutura de valores em ponto flutuante por um valor inteiro e um fator multiplicativo na base 10. Ele ainda insere na estrutura alguns bits contendo o sinal do valor e uns flags para casos especiais. Algo assim:

__attribute__((packed)) struct fp_s {
  unsigned int signalbit:1;  /* 0=positivo, 1=negativo */
  unsigned int factor:7;     /* fator multiplicativo de 10^(factor) */
  unsigned int flag:1;       /* 0=número válido, 1=número inválido */
  unsigned int number:23;    /* valor inteiro entre 0 e (2^24 - 1) */
};

A ideia é “deslocar a vírgula’ do valor inteiro para a esquerda ou direita do número inteiro, multiplicando-o pela potência de 10. O flag está ai para indicar se o “número” é infinito ou inválido. O que ele provavelmente não percebeu é que isso é semelhante a estrutura de um float, exceto pelos flags e que o expoente é computado através de potência de dois ao invés de dez. No caso do padrão IEEE 754 os valores zero, infinitos e “not a number” são obtidos com codificações especiais no container da variável. E, se a codificação decimal é tão importante assim, o padrão IEEE 754, na sua revisão de 1985, implementa os containers decimal32 e decimal64 e decimal128. Os tipos _Decimal32, _Decimal64 e _Decimal128 são suportados pelo GCC e pela especificação ISO C11, mas a biblioteca padrão não fornece suporte adequado para esses tipos!

Por que, num float, o expoente refere-se a uma potência de dois? Bem… dessa forma, ao invés de multiplicarmos ou dividirmos o valor inteiro por 10^e podemos, simplesmente, fazer shifts para esquerda ou direita do valor binário contido em number… E, para tornar as coisas ainda mais interessantes, o valor em number corresponde apenas à parte fracionária do valor, adicionando um 1.0 (binário) implícito e, depois, deslocando o ponto de acordo com o expoente… É por esse motivo que, num float, o valor 0.5, nos 4 bytes que compõem o valor armazenado, podem ser mapeados ao unsigned int 0x3f000000:

\displaystyle \begin{matrix}  \underbrace{0} & \underbrace{01111110} & \underbrace{0000\cdots000} \\  sinal & expoente & mantissa\,(23\,bits)  \end{matrix}

O exemplo abaixo mostra isso:

$ gcc -x c - <<@@@
#include <stdio.h>
float x = 0.5; 
void main(void)
{ printf("0x%08x\n", *(unsigned int *)&x); }
@@@
$ ./a.out
0x3f000000

O expoente 0x7E (126, veja a disposição dos bits acima) nos diz que o ponto binário tem que ser deslocado uma posição para a esquerda. Assim, 1.0 (2^0) torna-se 0.5 (2^{-1}).

O problema fundamental do ponto flutuante é que operações aritméticas feitas com eles nem sempre são exatas. Mas, temos os mesmos problemas com aritmética na base 10, se as fizermos manualmente. Como é fácil demonstrar, \frac{1}{3} não pode ser expresso exatamente sem usarmos a notação de fração. No entanto, a codificação binária dos float‘s e double‘s nos dão a possibilidade de armazenamento de gualquer valor inteiro que possa caber na mantissa mais 1 bit setado à esquerda de forma exata!

Pense bem… no caso do float temos 23 bits de precisão mais o bit implícito que compõe a unidade de valor. Isso nos permite o armazenamento exato de 24 bits. No caso do double, temos 53 bits, pelo mesmo motivo. E isso é feito automaticamente, através de instruções especiais do processador e das conversões feitas pelo próprio compilador.

Valores inteiros que cabem exatamente num float estão entre 0 e ±16777215 (24 bits). Para double‘s, entre 0 e ±9007199254740991 (53 bits, pouco mais de 9 quatrilhões). Para long long double temos a faixa toda dos 64 bits, de 0 até ±18446744073709551615 (pouco mais de 18.4 quintilhões). O desenvolvedor pode ficar tentado a usar o tipo long double todas as vezes, mas ele tem um sério problema. Extensões do processador, como SIMD, não suportam esse tipo. Ele é computado, necessariamente, usando instruções do co-processador matemático e gasta 10 bytes de espaço. Vejam a diferença de duas rotinas simples:

double f(double a, double b) { return a + b; }
long double g(long double a, long double b) { return a + b; }

No primeiro caso, apenas uma instrução é usada. No segundo, o compilador é obrigado a empilhar os valores, carregá-los para depois realizar a adição. São 3 instruções.

f:
addsd xmm0,xmm1
ret

g:
fld tbyte [rsp+8]
fld tbyte [rsp+24]
faddp st1, st
ret
Anúncios