O problema do ponto-flutuante (parte 3)

Dei um susto em vocês, né não? Que história é essa que ponto flutuante não atente a alguns axiomas da matemática? A lei da associatividade não funciona? A distributiva também não?! Pois é! É isso, meninos e meninas, quando se lida com ponto flutuante algumas leis universais da realidade vão por água abaixo!

A coisa piora quando falamos de erros relativos (“relativos” ao valor original!). Aquela coisa de que alguns números não podem ser representados com exatidão. É o exemplo de 0.1! Para demonstrar isso, eis um programinha:

/* f2u32.c */
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>

void showFloat(float x)
{
  int s, e, f;
  uint32_t u32 = *(uint32_t *)&x;

  /* Isola o sinal. */
  s = (u32 >> 31) & 1;

  /* Isola o expoente */
  e = (u32 >> 22) & 0xff;

  /* Isola o valor fracionário */
  f = u32 & 0x7fffff;

  printf("%f: 0x%08X {s:%d, e:0x%02X, f:0x%06X\n",
    x, u32, s, e, f);
}

int main(int argc, char **argv)
{
  float x;

  if (argc != 2)
  {
    fprintf(stderr, "Usage: f2u32 \n");
    return 1;
  }

  x = (float)atof(argv[1]);
  showFloat(x);

  return 0;
}

A rotina acima mostra a estrutura de um float que é passado na linha de comando. Se chamarmos o programa passando o valor 0.1, obtemos o seguinte:

$ ./f2u32 0.1 0.100000: 0x3DCCCCCD {s:0, e:0xF7, f:0x4CCCCD

Repare que a parte fracionária possui um período. O valor 0x4CCCCD, em binário, fica: 10011001100110011001101 (ignorando o bit 0 superior, já que só temos 23 bits!). Os últimos 3 bits (à direita) são arredondados para cima com base no padrão 1001 (1001₂ + 1₂ =1010₂).

Aumentar a precisão não adianta. Se usarmos o tipo double obteremos algo muito parecido (altere o código acima e tente!). O tipo long double tem expoente de 15 bits e parte fracionária de 64 bits. O valor equivalente de 0.1, em binário, no formato long double, é 0.1000000000000000055511151231257827. Esse restinho ai é graças à imprecisão de 2^-64 (o bit 0, mais à direita).

Só para ilustrar, o gráfico abaixo mostra o nível de precisão dos números representados num float e o conjunto dos números reais:

Pontos mais altos são mais exatos que pontos mais baixos!

Ok. O gráfico não tem escala e eu mesmo acho que não está correto, mas foi o único gráfico que achei para demonstrar que apenas um subconjunto estreito dos números reais pode ser representado com precisão considerável em ponto flutuante. De forma geral, qualquer valor real que possa ser representado pela somatório de fatores 2^-n, onde n seja menor ou igual ao número de bits de f, pode ser representado com exatidão.

A figura abaixo mostra a estrutura de um float em detalhes:

O valor do expoente (e) tem alguns significados:

Se e for igual a zero e a fração (f) for zero, então temos uma anomalia interessante: Dependendo do sinal (s) podemos ter +0 ou -0;

Se e for igual a 255 e a fração (f) for zero, dependendo do sinal, temos +infinito e –infinito. Mas se f for diferente de zero, temos +NaN ou –NaN. A sigla NaN significa “Not a Number“.

Valores válidos para e, com um valor diferente de zero para f estão compreendidos na faixa de 1 até 256. O valor zero para e com f diferente de zero também tem um significado especial. São valores chamados de subnormal. Não os discutirei aqui.

Mas, pera ai… tem alguma coisa errada! Se o número 0.1 não pode ser representado em ponto flutuante, como é que printf() consegue mostrá-lo corretamente?! Mágica?! Acontece que as rotinas que usam ponto flutuante trabalham com uma margem de erro maior do que o erro do tipo float (ou double, ou long double). Para demonstrar isso, veja o seguinte fragmento de código:

float x;
uint32_t *p = (uint_t *)&x;

for (*p = 0x3DCCCC89; *p < 0x3DCCCD11; *p++)
  printf("%f: 0x%08X\n", x, *p);

Você verá que se fizermos *p variar de 0x3DCCCC8A até 0x3DCCCD10, obteremos o mesmo valor 0.10000, no printf(). Se você estipular a precisão no parâmetro %f, diminuirá o erro interpretado pela função.

Quero que você entenda que a precisão com que printf() mostra o valor decimal nada tem haver com a precisão de uma variável float. A função printf() sabe como decodificar um float e formatá-la de forma a parecer que tem a precisão que você desejar. Lembre-se a formatação ‘f’ segue o padrão:

%[width][.precision]f

Onde width é quantidade de algarismos e precision é a quantidade de “digitos decimais”, após o ponto decimal. No exemplo acima, se aumentarmos a precisão para 15, por exemplo, e inicializarmos o conteúdo do ponteiro p para 0x3DCCCCCC e 0x3DCCCCCD, obteremos 0.0999999940 e 1.0000000015, respectivamente.

O que meio que reforça o meu ponto: 0.1 não pode ser representado em ponto flutuante! Mas não impede de printf() format o valor com precisão!

Repare, no entanto, que embora printf() mostre o mesmo resultado para a precisão default, os valores são diferentes entre si. A minha recomendação padrão para valores em ponto flutuante segue a piada de Grouxo Marx: “Você vai confiar em mim ou nos seus próprios olhos?” – Neste caso, confiar nos seus olhos pode gearar alguns problemas! (Confiar em mim também pode! Mas, hey, você sempre pode verificar as coisas por si mesmo, não pode?).

Voltando ao erro relativo. Como vímos, podemos ter um erro de 2^-23 (correspondendo ao último bit, mais à direita, da fração “f”) no valor final atribuído a uma variável float. Esse erro é agravado quando realizamos operações de adição e subtração. Nessas operações o erro é acumulado até produzir um erro maior que o esperado. O mesmo acontece com multiplicações e divisões, mas numa escala menor, isto é, multiplicações e divisões propagam menos erros que adições e subtrações… Não vou demonstrar isso aqui porque tomaria muito espaço no texto. Para uma referência (que não é definitiva) sobre esses erros, consulte o volume 2 da maravilhosa série de Dr. Donald E. Knuth, The Art of Computer Programming… Lá ele não trata do padrão IEEE 754, mas é um embasamento matemático muito interessante sobre aritmética de ponto flutuante.

Graças a esse pequenino erro temos o problema das leis associativa e distributiva indo pro saco… Felizmente a lei comutativa não sofre esse problema:

a + b = b + a
a * b = b * a

Outras leis também são garantidas:

a – b = a + -b = -b + a

Leia o primeiro capítulo do livro do Dr. Knuth e veja o que é seguro e o que não é (de novo, longe de ser uma discussão definitiva sobre o assunto, ok?).

Embora esses axiomas, ou essas leis, sejam válidas no uso de ponto flutuante, deve-se levar em consideração a propagação dos erros. É por causa deles que, às vezes, quando você adiciona R$ 0.99 a R$ 4.01, por exemplo, o resultado é R$ 5.01.

Você deve estar pensando: “E como é que os bancos de dados, como o Oracle, por exemplo, lidam com ‘números quebrados’?”. A declaração de valores decimais é feita assim:

price decimal(10,2) not null

O princípio aqui é o mesmo que o usado em printf(). O valor tem 10 digitos (contando com os “depois da vírgula”), sendo que a precisão é de 2 digitos “depois da vírgula”. Os DBMSs costumam armazenar esses valores em decimal codificado em binário (BCD). Em BCD cada conjunto de 4 bits armazena um algarismo decimal (entre 0000₂ e 1001₂). O DBMS então pode controlar como ele trata o erro residual das multiplicações e divisões (somas e subtrações não geram erros, neste caso!). Por isso, se você obtém o valor de cum campo decimal, do banco de dados, e o armazena numa variável do tipo float, double ou long double com o intuito de fazer cálculos, pode ser que obtenha erros que o banco de dados não daria. A especificação de C++0x (TR1) prevê o a definição de tipos decimais, justamente para ocupar essa lacuna nas especificações de C e C++. O “dó sustenido”, C#, e o idioma Pascal usado pelo Delphi já possuem um tipo decimal (no Delphi, chama-se Currency)… Só que ambos são implementados de forma diferente. Eles são inteiros! Só que o valor que representam é dividido por 100 para representar 2 casas decimais no valor (isso é um exemplo – não me lembro se Currency é dividido por 100 ou por 10000, se são 2 ou 4 casas decimais.

Aritmética usando BCD é assunto para outro post. A especificação C++0x (N3126) também é assunto para outro post.