Re: Floating point problem on glibc-2.1.1

[Tung-Han Hsieh <thhsieh@linux.org.tw>]

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