Floating-point directed roundings and optimization

Question

There is the following code in which the same expression is evaluated in different rounding modes:

#include <iostream>
#include <fenv.h>
#pragma STDC FENV_ACCESS ON
#define SIZE 8

double foo(double * a, double * b){
    double sum = 0.0;
    for(unsigned int i = 0; i < SIZE; i++) {
        sum+= b[i] / a[i];
    }
    return sum;
}

int main() {
    double a[]={127, 131, 137, 139, 149, 151, 157, 163};
    double b[SIZE];

    for(unsigned int i = 0; i < SIZE; i++){
        b[i] = i+1;
    }

    printf("to nearest:   %.18f \n", foo(a, b));

    fesetround(FE_TOWARDZERO);
    printf("toward zero:  %.18f \n", foo(a, b));

    fesetround(FE_UPWARD);
    printf("to +infinity: %.18f \n", foo(a, b));

    fesetround(FE_DOWNWARD);
    printf("to -infinity: %.18f \n", foo(a, b));

    return 0;
}

When it is compiled with g++ with the -O0 option, the output is as follows:

to nearest:   0.240773868136782450
toward zero:  0.240773868136782420
to +infinity: 0.240773868136782560
to -infinity: 0.240773868136782420

But when compiling it with the -O3 option, we have:

to nearest:   0.240773868136782480
toward zero:  0.240773868136782480
to +infinity: 0.240773868136782480
to -infinity: 0.240773868136782480

Compiler: g++ (MinGW.org GCC-6.3.0-1) 6.3.0

Why don't the rounding modes change? How to fix it?

(if fesetround is called at each iteration of the for loop (inside the foo function), then the results are correct with any compilation flags.)

UPD: I think the problem is that the compiler computes the value of fesetround at compile type as @haneefmubarak has pointed out in https://stackoverflow.com/a/26319847/2810512. The question is how to prevent it. (only for one command, fesetround, not for the whole function).

I wrote wrappers for rounding routines with __attribute__ ((noinline)) and call them in the main function:

void __attribute__ ((noinline)) rounddown(){
    fesetround(FE_DOWNWARD);
}

void __attribute__ ((noinline)) roundup(){
    fesetround(FE_UPWARD);
}

int main() {
    ...

    roundup();
    printf("to +infinity: %.18f \n", foo(a, b));

    rounddown();
    printf("to -infinity: %.18f \n", foo(a, b));
    ...
}

But it does not work. Any ideas?

UPD2: A more clear example:

Correct rounding (-O0)

Failed rounding (-03)

It is easy to see that the exact result:

2/3 + 2/5 + 4/7 + 4/11 = 2.0017316017316017316...

If you store `foo(a,b)` in a variable `c` and use that variable in all your `printf` statements, do you have the same problem?. — kvantour, Oct 23 '18 at 14:03
Does the version of GCC you are using support `#pragma STDC FENV_ACCESS ON`? — Eric Postpischil, Oct 23 '18 at 14:40
@EricPostpischil there is a warning (g++ from MinGW) `pragma STDC FENV_ACCESS ON is not supported, ignoring pragma [-Wunknown-pragmas]` But it works with the `O0` flag. — Konstantin Isupov, Oct 24 '18 at 07:02
Could you provide us with some more information. OS, g++ version, compilation line. — kvantour, Oct 24 '18 at 07:17

score 5 · Accepted Answer · answered Oct 24 '18 at 10:25

Per a comment from the question author, the compiler they are using does not support #pragma STDC FENV_ACCESS ON and prints a warning saying so.

The code likely “works” in the unoptimized version because fesetround does change the rounding mode in the hardware and the compiler emits straightforward code doing operations in the nominal order represented by the source code.

Reasons the optimized code does not work may include:

The compiler performs some arithmetic at compile time, ignoring the fesetround calls.
During optimization, the compiler reorders operations, possibly performing arithmetic operations and fesetround calls in a different order than shown by the source code. The fesetround calls may even be removed completely.

There may be no fix for this in C. If the compiler does not support accessing the floating-point environment, there may be no way to force it to generate the necessary code. Declaring some objects volatile may force some operations to be performed at execution time and in a desired order, but a compiler might still reorder fesetround with respect to these operations, depending on what information about fesetround is built into it.

It may be necessary to use assembly language to perform floating-point arithmetic with desired rounding modes.

Floating-point directed roundings and optimization

1 Answers1