Os perigos do ponto flutuante para o Cálculo

Já escrevi algumas vezes aqui sobre ponto-flutuante… Aqui vai o alerta, de novo: Cálculo com diferenciais, em ponto-flutuante, podem causar erros. Eis um teste simples: Suponha que você tenha uma função y=x^2 e execute o seguinte teste:

/* square.c */
#include <stdio.h>

#ifdef FLOAT
  #define REAL float
  #define VALUE 2.0f
  #define MIN 5e-6f
  #define PREC "8"
#else
  #define REAL double
  #define VALUE 2.0
  #define MIN 5e-14
  #define PREC "16"
#endif

#define F(x) ((x) * (x))

void main(void)
{
  REAL x, dx, xpdx;
  REAL y, dy, s1r, s2r;

  x = VALUE;
  dx = MIN;

  // Calcula cada um dos componentes separadamente.
  y = F(x);
  dy = F(dx);
  s1r = y + dy;

  // Calcula de uma vez só.
  xpdx = x + dx;
  s2r = F(xpdx);

  printf("x = %." PREC "f\n"
         "dx = %." PREC "f\n"
         "x + dx = %." PREC "f\n"
         "f(x) = %." PREC "f\n"
         "f(dx) = %." PREC "f,\n"
         "f(x) + f(dx) = %." PREC "f\n"
         "f(x + dx) = %." PREC "f\n",
         x, dx, xpdx,
         y, dy, s1r, s2r);
}

No programinha acima fazemos x^2+dx^2 e (x + dx)^2 e mostramos os dois resultados. O valor dx é escolhido com sendo próximo ao limite da precisão do tipo escolhido (float ou double) intencionalmente, já que na teoria dx deveria aproximar-se tão próximo quando possível de zero. Em ambos as operações deveriamos obter o mesmo resultado, já que os dois métodos são equivalentes, mas isso não acontece:

$ gcc -O0 -DFLOAT -o squaref square.c
$ ./squaref
x = 2.00000000
dx = 0.00000500
x + dx = 2.00000501
f(x) = 4.00000000
f(dx) = 0.00000000,
f(x) + f(dx) = 4.00000000
f(x + dx) = 4.00002003

Pegue uma calculadora e você verá que f(x)+f(dx) é o resultado que mais se aproxima do valor verdadeiro:

\displaystyle \left(2+\frac{5}{10^6}\right)^2=2^2+\left(\frac{5}{10^6}\right)^2=4+\frac{25}{10^{12}}=4.000000000025

O que está acontecendo, então, quando fazemos f(x+dx)? E, por que f(dx) = 0.0?

Em meu livro (aqui) explico que os valores em ponto flutuante são codificações em “notação científica”, mas em base 2 ao invés da badalada base 10. Isso quer dizer que todo número em ponto flutuante tem uma “mantissa” e um expoente com quantidade de bits limitados.

Acontece que adições (e subtrações) de valores em notação científica têm que levar em conta o expoente aplicado à base… Quando você faz 2 + 5 \cdot 10^{-1}, por exemplo, está fazendo 2 \cdot 10^0 + 0.5 \cdot 10^0. Como os expoentes são iguais você pode realizar a adição e obter 2.5 \cdot 10^0. Repare que o valor menor é deslocado para a esquerda. E, como temos um número limitado de bits de precisão na mantissa, quanto mais deslocado para a esquerda, maior é a chance desses bits, mais à esquerda, sumirem! Isso é chamado erro de truncamento.

No caso do tipo float a precisão da mantissa é de 24 bits (23 bits mais o bit inicial ‘1’, em valores “normalizados”), o que nos dá, em decimal, uma precisão de, no máximo, 7 algarismos (\log 2^{24} = 7). Quando fazemos dx^2, perdemos todos os bits de dx na multiplicação porque o expoente -6 torna-se -12 (e, nessa conta,  só podemos chegar a -6 porque temos apenas 7 alguarismos possíveis em um valor!). Por isso f(dx) resulta em 0.0. No outro caso, f(x+dx), a adição não causa perdas significativas porque os dois valores cabem em 24 bits. Mas teremos truncamento na multiplicação… Os processadores Intel fazem suas operações aritméticas em ponto flutuante usando uma precisão maior que o tipo definido no seu código (provavelmente usando o tipo long double, com mantissa de 64 bits), mas o resultado tem que ser “truncado” de volta para a precisão desejada, “comendo” alguns bits no processo…

Obviamente que podemos aumentar a precisão do resultado aumentando a quantidade de bits da mantissa, de 23 para 52 bits (ao invés de usar  float, usar double). Isso nos dará uma margem de trabalho de 15 algrarismos decimais, mas o truncamento de bits residuais continuará acontecendo. Subtituindo o tipo float por double (compilando sem a definição FLOAT) na rotina lá em cima (e também substituindo dx por 5e-14), teremos este resultado:

$ gcc -O0 -o square square.c
$ ./square
x = 2.0000000000000000
dx = 0.0000000000000500
x + dx = 2.0000000000000502
f(x) = 4.0000000000000000
f(dx) = 0.0000000000000000,
f(x) + f(dx) = 4.0000000000000000
f(x + dx) = 4.0000000000002007

Embora pareça que f(x+dx) dê um resultado mais aproximado, ele continua sofrendo do mesmo problema… ele “comeu” quase 2.3 \cdot 10^{-11} do valor correto! A aparência de melhor aproximação é enganosa, neste caso, porque esse valor “depois da vírgula” é devido a erro de truncamento, do mesmo jeito que no float. A operação f(x)+f(dx), neste caso específico, continua sendo mais precisa, mesmo que tenha um erro maior!

Erros de truncamento não são o único problema… Quando uma operação em ponto flutuante é feita, os bits inferiores à precisão desejada são usados para “arredondar” o valor final… É como se, em decimal, pudéssemos usar apenas valores inteiros e uma conta dê o resultado 1,7… Esse valor 1.7 pode ser arredondado para 2.0, causando a perda do algarismo depois da vírgula… No caso dos floats o arredondamento pode ser feito (pelo processador, automaticamente) adicionando 5 \cdot 10^{-7} à mantissa. No caso de double, 5 \cdot 10^{-15}.

Meu ponto aqui, e sempre que falo sobre ponto flutuante, é o seguinte: Não adianta nada você achar que só porque tem a “precisão” de algarismos “depois da vírgula” terá, exatidão nos cálculos… Em alguns casos o valor obtido estará longe de ser o valor real. É absolutamente necessário levar em conta o funcionamento das operações em ponto flutuante e, além disso, aquilo que os matemáticos chamam de estabilidade no cálculo. Especialmente se os dois argumentos tiverem grandezas muito discrepantes!

Anúncios

Deixe um comentário

Faça o login usando um destes métodos para comentar:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair / Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair / Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair / Alterar )

Foto do Google+

Você está comentando utilizando sua conta Google+. Sair / Alterar )

Conectando a %s