5 replies to this topic

[monjardin]

Can anyone explain why this code:

#include <stdio.h>


int main(int argc, char *argv[])

{

        double a = -0.6;

        double b = 2025.0;

        double c = a * b;


        printf("a * b = %f\n", a * b);

        printf("a * b = (int) %d\n", (int) (a * b));

        printf("a * b = (float) %f\n", (float) (a * b));

        printf("a * b = (int) (float) %d\n", (int) (float) (a * b));

        printf("a * b = (int) %d\n", (int) c);


        return 0;

}

produces this result:

a * b = -1215.000000

a * b = (int) -1214

a * b = (float) -1215.000000

a * b = (int) (float) -1215

a * b = (int) -1215

Of particular interest is the second integer result of 1214. The results were produced with gcc v3.4.3.

[Reedbeta]

Well, I'd guess that in the second case, a*b is being computed with 80 bit internal FPU precision (assuming you're using x86) and converted to integer from there, while #4 and #5 involve the product being rounded down to 64 bits or fewer before being converted to integer.

Conversion to integer is done by truncation of the fractional part, so perhaps the 80 bit result contains a number slightly greater than -1215, i.e. -1214.9999.... which is truncated to -1214 when converted to integer directly, but rounded to exactly -1215 when converted to double or float.

Just my guess

[monjardin]

That sounds good to me. It is an Intel chip.

A coworker of mine was tracking down a bug and narrowed it down to this rounding error. I'm used to being able to explain these events, but this one had me puzzled. It only occurred when the value was never stored in a variable and cast directly to an integer. I even had him print the 64-bit hexadecimal value so I could convert it by hand. ;)

Thanks.

[Reedbeta]

When dealing with floats I frequently find these three pages quite useful.

By the way, in case you didn't know, a nice trick to get rounding to nearest integer (instead of trucation) is add 0.5 if it's positive, or subtract if it's negative, before casting.

[Nick]

You can adjust the internal precision with the x87 FPU control word. In Visual C++ with inline assembly you can do it like this:

void setSinglePrecision()   // 32-bit

{

    unsigned short cw;


    __asm

    {

        fstcw cw

        and cw, 0xFCFF

        fldcw, cw

    }

}


void setDoublePrecision()   // 64-bit

{

    unsigned short cw;


    __asm

    {

        fstcw cw

        and cw, 0xFEFF

        or cw, 0x0200

        fldcw, cw

    }

}


void setDoubleExtendedPrecision()   // 80-bit

{

    unsigned short cw;


    __asm

    {

        fstcw cw

        and cw, 0xFDFF

        or cw, 0x0100

        fldcw, cw

    }

}

I'm curious whether that 'fixes' your issue though...

[monjardin]

Reedbeta said:

When dealing with floats I frequently find these three pages quite useful. :)

I'm very familiar with those! I used the second one to very that the answer is not off by some small value with a 64-bit double. Then I came here. ;)

Nick said:

fstcw cw
and cw, 0xFEFF
or cw, 0x0200
fldcw, cw

I'm out of the office today, but I'll email in this snippet and see what happens.