Floating point problem on glibc-2.1.1
Tung-Han Hsieh
thhsieh@linux.org.tw
Tue May 4 01:11:00 GMT 1999
On Mon, May 03, 1999 at 08:38:51AM -0700, Ulrich Drepper wrote:
> The problem is in you code. What do expect other than an unreliable
> result? All implementations have some error and if it is (in this
> case) the case that the error struck at the wrong point you get
> surprising results. If you would use
>
> printf("%ld\n", round(log(8.0)/lrint(2.0)));
>
> instead you'd get the expected result.
Hello,
I try the code you mentioned, but I got the result "524288" in glibc-2.1.1,
too far from the result I expected.
In fact, I heard that this problem was reported from the Fermi Lab, and
they used the code like the following for test:
#include <stdio.h>
#include <math.h>
main()
{
int j, j1;
double dj, d=8.0;
dj = (1.0 + log(d) / log(2.0));
j = (int) dj;
j1 = (int) (1.0 + log(d) / log(2.0));
printf("j=%d, j1=%d, dj=%f\n", j, j1, dj);
}
In glibc-2.1.1, the result is
j=4, j1=3, dj=4.000000
but we know that the value of "j1" is not correct. So you mean this is
the unreliable result?
Here is another test program written by <dinon.bbs@bbs.ee.ncu.edu.tw>:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef unsigned long long verylong;
// sign: 1, exponent: 11, mantissa: 52
void hex(double d)
{
int i;
verylong val;
unsigned int part[2];
memcpy(&val, &d, 8);
memcpy(part, &d, 8);
printf("val = %lf hex = %016llX ", d, val);
printf("s = %d e = %d ",
(part[1] >> 31),
(part[1] & 0x7FF00000) >> 20);
// extract the most sigifinicant 20 bits of mantissa
part[1] &= 0xFFFFF;
memcpy(&val, part, 8);
printf("m = %016llX\n", val);
}
int main()
{
int j, j1;
double dj, d = 8.0;
dj = (1.0 + log(d) / log(2.0));
j = (int) dj;
j1 = (int) (1.0 + log(d) / log(2.0));
printf("j=%d j1=%d dj=%f\n", j, j1, dj);
hex(4.0);
hex(j1);
return 0;
}
When I compile it under glibc-2.1.1, I got the result:
j=4 j1=3 dj=4.000000
val = 4.000000 hex = 4010000000000000 s = 0 e = 1025 m = 0000000000000000
val = 3.000000 hex = 4008000000000000 s = 0 e = 1024 m = 0008000000000000
Does glibc treat this phenomena as "problem" or "feature"?
(PS. Sorry, I might forget to mention. I use Pentium-Pro 200 and Pentium-133
(no MMX) for test.)
T.H.Hsieh
More information about the Libc-alpha
mailing list