unexpected rounding behaviour of mpf_add(). bug?

[Frederik Schaffalitzky]

I use the mpf_t data type with a fixed precision (256 bits everywhere).
While I can estimate the "machine" epsilon as the largest power of 1/2
such that
  1 + epsilon == 1

holds exactly (as deemed by mpf_cmp()), for many numbers `a' there is no
power of 1/2 which has the property that
  a + epsilon == a

holds exactly. E.g. with 64 bits of precision, and a trite choice of `a',
the result of "b = a + eps" is:

a: 0.14142135623730950488e1
b: 0.141421356237309504876e1
eps: 0.120061587716658067602e-1671

a: 0.14142135623730950488e1
b: 0.141421356237309504876e1
eps: 0.600307938583290338012e-1672

a: 0.14142135623730950488e1
b: 0.141421356237309504876e1
eps: 0.300153969291645169006e-1672

etc etc, regardless of how small I make `eps'. I.e. adding a tiny number
to `a' seems to introduce significant (binary) digits a long way away from
the significant (binary) digits of that tiny number.

I am not sure if this is strictly speaking a bug as GMP makes few
guarantees about rounding behaviour, but on the other hand it is hard to
imagine an addition algorithm which produces such results and which isn't
rather careless about trailing digits. Is there some way I can make mpf_t
behave in the way that I desire?

[Whatever the merit of this practice, I have a numerical algorithm which
assumes this sort of behaviour for its termination criterion.]

Thank you,

Frederik

--------------------------------------------------------------------

My GMP is version 3.1.1; I built it from the tarfile using
  ./configure --cache-file=/dev/null --prefix=`cd ../install; pwd`


$ gcc -v
Reading specs from
/data/az0/fsm/gcc-2.95.2/install/lib/gcc-lib/i686-pc-linux-gnu/2.95.2/specs
gcc version 2.95.2 19991024 (release)


$ uname -a
Linux shadow.robots.ox.ac.uk 2.2.17 #4 SMP Tue Feb 20 13:28:13 GMT 2001
i686 unknown


$ ./config.guess 
pentium3-pc-linux-gnu


How to reproduce:
Compile and link the test program against libgmp. Run with no arguments.
Symptom: the program does not terminate.

-----------------------------------------------------

#include <stdio.h>
#include <gmp.h>

int main()
{
  int N = 64;
  mpf_t a, b, eps, factor;
  
  /* a = sqrt(2) */
  mpf_init2(a, N);
  mpf_set_d(a, 2.0);
  mpf_sqrt(a, a);
  
  /* ... */
  mpf_init2(b, N);
  
  /* eps = 1 */
  mpf_init2(eps, N);
  mpf_set_d(eps, 1.0);
  
  /* factor = 2 */
  mpf_init2(factor, N);
  mpf_set_d(factor, 2.0);
  
  while (1) {
    /* eps /= factor */
    mpf_div(eps, eps, factor);
    
    /* b = a + eps */
    mpf_add(b, a, eps);
    
    printf("a: "); mpf_dump(a);
    printf("b: "); mpf_dump(b);
    printf("eps: "); mpf_dump(eps);
    printf("\n");
    
    if (mpf_cmp(a, b) == 0)
      break;
  }
  
  return 0;
}