a flirt with epsilon

[Francois Grieu]

When I run the following code onsome (obscure) platforms,
I get strange results.


#include <stdio.h>
#include <float.h>
#include <math.h>

/* attempt to find something like DBL_EPSILON at run time */
double dbl_epsilon(void)
{
double x = 1, y;
do
{
y = x;
x *= 0.5;
}
while (x+1!=1);
return y;
}

int main(void)
{
#define log14_13 0.074107972153721878469097423336037692907
printf("DBL_EPSILON %g\n",DBL_EPSILON);
printf("dbl_epsilon %g\n",dbl_epsilon());
printf("log rel err at 14/13 %g\n",log(14./13)/log14_13-1);
printf("log1p rel err at 1/13 %g\n",log1p(1./13)/log14_13-1);
return 0;
}


On one platform I get:
DBL_EPSILON 2.22045e-16
dbl_epsilon 1.0842e-19 <- way better than DBL_EPSILON
log rel err at 14/13 5.42101e-19
log1p rel err at 1/13 -1.0842e-19


DBL_EPSILON 2.22045e-16
dbl_epsilon 2.22045e-16
log rel err at 14/13 -5.95187e-16
log1p rel err at 1/13 -3.33936e-17

How can the later result be much better than DBL_EPSILON,
by over 3 bits ?

TIA,

François Grieu