Aritmética com ponto-flutuante arbitrária

1. Motivação
    A aritmética de ponto-flutuante utilizada nos computadores é susceptível a alguns erros relacionados  com a forma que os números reais são representados (base binária) e também relacionado com o tamanho da representação (32 bits ou 64 bits). Esses dois fatores limitantes frequentemente geram erros de precisão nos cálculos realizados com os números com ponto flutuante.


    Como na base decimal, certos números reais tem uma representação infinita de dígitos. Mesmos alguns números com a representação na base decimal finita, podem ter uma quantidade infinita de dígitos na base binária. Como, por exemplo,  o número decimal 0.1 tem representação binária 0.0(0011). Os números entre parênteses repete-se indefinidamente. Por outro lado, o número decimal 0.125 tem representação de ponto-flutuante finita na base binária 0.001.


  Um número ponto flutuante é armazenado da seguinte maneira (+|-) s*2^e. O s representa o sginificante. O significante pode ter 24 ou 53 dígitos binários.


    Por exemplo, o número 0.1 é representado no computador da seguinte maneira:


 1.1001100110011001100110011001100110011001100110011010*2^-4


    Essa representação ocasionará erros perceptíveis em cálculos simples:


    Considere o seguinte programa em C:

 double x = 1.0;
 for(i=0;i<100;i++){
 x = x - 0.01;
 } 
 printf("%.15lf\n",x);
 

   A saída desse programa será:   
  -0.000000000000001




    A idéia apresentada em [1] consiste em calcular o erro de precisão associado uma operação armazená-la e propagá-la com o resultado do cálculo. Dessa maneira, o resultado de uma operação deixa de ser apenas um ponto-flutuante e passa a ser um expansão e = en+en-1+...+e2+e1, onde cada componente ei é um número ponto-flutuante.


    A operação de adição de dois números de ponto flutuante fica da seguinte maneira:

#define Fast_Two_Sum_Tail(a, b, x, y) \
 bvirt = x - a; \
 y = b - bvirt
 #define Fast_Two_Sum(a, b, x, y) \
 x = (REAL) (a + b); \
 Fast_Two_Sum_Tail(a, b, x, y)

  
 A expansão e gerada pela adição de a e b é x+y.
  
  Figura extraída qualificação de doutorado Markos Oliveira Freitas


 A soma de expansão com um número com ponto-flutuante fica da seguinte maneira:



    Figura extraída qualificação de doutorado Markos Oliveira Freitas


     Utilizando essa idéia, o algoritmo fica da seguinte maneira:   

  int grow_expansion(elen, e, b, h) /* e and h can be the same. */
  int elen; //tamanho da expansão
  REAL *e;  //expansão
  REAL b;   //ponto-flutuante adicionado
  REAL *h;  //expansão de saída

  e[0] = 1.0;
  h[0] = 0.0;
  for(i=0;i<100;i++){
   elen = grow_expansion(elen,e,-0.01,h);
   expansioncopy(elen,e,h);
  }
  printf("%.15lf\n", estimate(elen, e) );

O resultado do programa será  

0.000000000000000

O código C com as operações pode ser encontada aqui.




Referências
[1] J. R. Shewchuk. Adaptive precision floating-point arithmetic and fast robust   geometric predicates. Discrete & Computational Geometry, 18(3).