A flirt with epsilon

[Francois Grieu]

When I run the following code ons ome (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


On another platform I get:
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 nearly 3 bits ?

TIA,

François Grieu

 

[Martin Dickopp]

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

I can reproduce this even on a mainstream platform (x86, gcc 3.3).

#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


On another platform I get:
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 nearly 3 bits ?

I'm not sure your function to calculate DBL_EPSILON is correct. In
particular, the compiler can optimize the expression `x + 1 != 1'.

On my platform, if I replace said expression with `!cmp (x + 1.0, 1.0)'
and define

int cmp (double x, double y) { return x == y; }

in a different translation unit (to defeat any optimization attempt),
`dbl_epsilon ()' yields 2.22045e-16.

Martin

 

[CBFalconer]

When I run the following code ons ome (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

On another platform I get:
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 nearly 3 bits ?

You will probably find that your loop is using something like 80
bit floats, and leaving all results in the FP processor.
Investigate your compiler flags.