c++ - Why does returning a floating-point value change its value?

The following code raises the assert on Red Hat 5.4 32 bits but works on Red Hat 5.4 64 bits (or CentOS).

On 32 bits, I must put the return value of millis2seconds in a variable, otherwise the assert is raised, showing that the value of the double returned from the function is different from the one that was passed to it.

If you comment the "#define BUG" line, it works.

Thanks to @R, passing the -msse2 -mfpmath options to the compiler make both variants of the millis2seconds function work.

/*
 * TestDouble.cpp
 */

#include 
#include 

#include 

static double millis2seconds(int millis) {
#define BUG
#ifdef BUG
    // following is not working on 32 bits architectures for any values of millis
    // on 64 bits architecture, it works
    return (double)(millis) / 1000.0;
#else
    //  on 32 bits architectures, we must do the operation in 2 steps ?!? ...

    // 1- compute a result in a local variable, and 2- return the local variable
    // why? somebody can explains?
    double result = (double)(millis) / 1000.0;
    return result;
#endif
}

static void testMillis2seconds() {
    int millis = 10;
    double seconds = millis2seconds(millis);


    printf("millis                  : %d\n", millis);
    printf("seconds                 : %f\n", seconds);
    printf("millis2seconds(millis)  : %f\n", millis2seconds(millis));
    printf("seconds <  millis2seconds(millis)  : %d\n", seconds < millis2seconds(millis));
    printf("seconds >  millis2seconds(millis)  : %d\n", seconds > millis2seconds(millis));
    printf("seconds == millis2seconds(millis)  : %d\n", seconds == millis2seconds(millis));

    assert(seconds == millis2seconds(millis));
}


extern int main(int argc, char **argv) {
    testMillis2seconds();
}

Answer

With the cdecl calling convention, which is used on Linux x86 systems, a double is returned from a function using the st0 x87 register. All x87 registers are 80-bit precision. With this code:

static double millis2seconds(int millis) {
    return (double)(millis) / 1000.0;

};

The compiler calculates the division using 80-bit precision. When gcc is using the GNU dialect of the standard (which it does by default), it leaves the result in the st0 register, so the full precision is returned back to the caller. The end of the assembly code looks like this:

fdivrp  %st, %st(1)  # Divide st0 by st1 and store the result in st0
leave
ret                  # Return

With this code,

static double millis2seconds(int millis) {
    double result = (double)(millis) / 1000.0;
    return result;
}

the result is stored into a 64-bit memory location, which loses some precision. The 64-bit value is loaded back into the 80-bit st0 register before returning, but the damage is already done:

fdivrp  %st, %st(1)   # Divide st0 by st1 and store the result in st0
fstpl   -8(%ebp)      # Store st0 onto the stack
fldl    -8(%ebp)      # Load st0 back from the stack
leave
ret                   # Return

In your main, the first result is stored in a 64-bit memory location, so the extra precision is lost either way:

double seconds = millis2seconds(millis);

but in the second call, the return value is used directly, so the compiler can keep it in a register:

assert(seconds == millis2seconds(millis));

When using the first version of millis2seconds, you end up comparing the value that has been truncated to 64-bit precision to the value with full 80-bit precision, so there is a difference.

On x86-64, calculations are done using SSE registers, which are only 64-bit, so this issue doesn't come up.

Also, if you use -std=c99 so that you don't get the GNU dialect, the calculated values are stored in memory and re-loaded into the register before returning so as to be standard-conforming.