Machine Epsilon (Dr. Ing. Carlos L?pez)

[Jim White]

If we apply the standard definition of machine epsilon
to a "virtual machine" in the form of an arbitrary
precision arithmetic software system, then we should
theoretically have an arbitrarily selectable espilon,
and indeed this will vary from variable to variable
according to the sizes allocated to each.

This should be entirely predictable, as the GMP and
MPFR floating-point data types are fixed-precision
(this can be arbitrarily selected of course).

If we apply the standard test, "what's the maximum EPS
for which  (1 +  1 / 2^EPS) > 1  returns true", then
we would expect the epsilon value to be fixed by the
underlying arithmetic base.

For GMP's mpf_t type, on a 32-bit platform, if I have
a float variable with K limbs, then EPS is 32K.  It's
not 32K-1 because GMP allocates an extra limb as a
"guard digit" for mpf_t's.

If I have a 3-limb, 96-bit mantissa, then I can add 1
+ 1/2^95, and have a result exactly representable
within my 96-bit mantissa.

But I can also add 1 + 1/2^96 and still get a result
distinguishable from 1.0, but not 1 + 1/2^97.

So my variable's epsilon is 96.

Pseudo-code for verification:

   int eps = 0;
   mpf_t x, y;
   mpf_init2 (x, 96);
   mpf_init2 (y, 96);

   mpf_set_ui (y, 1);
   do {
      mpf_div_ui (y, 2);
      mpf_add_ui (x, y, 1);
      if (mpf_cmp_ui(x,1) == 0) break;
      eps++;
      }


MPFR uses GMP's integer arithmetic services to
implement its own floating-point type, mpfr_t.  The
main difference is that MPFR allows mantissa
precisions to be specified in bits, not whole limbs,
and provides a choice of rounding schemes.

Changing the code above to mpfr, with "round to
nearest", the epsilon value for a mantissa of 96 bits
is returned as 95, which is exactly what we would
expect. 

The internal floating-point representations of both
systems are well defined in the respective manuals. 


Cheers

Jim White
ANU Canberra