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]

Francois Grieu <fg****@francenet.fr> writes:

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
--
,--. Martin Dickopp, Dresden, Germany ,= ,-_-. =.
/ ,- ) http://www.zero-based.org/ ((_/)o o(\_))
\ `-' `-'(. .)`-'
`-. Debian, a variant of the GNU operating system. \_/

[CBFalconer]

Francois Grieu wrote:

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.

--
"I'm a war president. I make decisions here in the Oval Office
in foreign policy matters with war on my mind." - Bush.
"Churchill and Bush can both be considered wartime leaders, just
as Secretariat and Mr Ed were both horses." - James Rhodes.