How to prevent C99 floating point code from changing results with optimization level

Question

I made this code below to test a system's compliance with C99's IEEE 754. The problem is the results change with different optimization levels. If I stick with optimization level zero, most of the tests pass. But if I increase the optimization level to -O3, then most of the tests fail. I have tried -ffloat-store. It didn't change any results.

I have tried this program on Mac OS 10.4.11 PowerPC using GCC 5.2.0 and Mac OS 10.6.8 using GCC 4.2.1 and GCC 4.9.2. The results are the same. Most of the tests pass only when using optimization level 0 (-O0). I am thinking this is a bug with GCC, but before I report it, I would like to hear others' opinions on the code.

This is how I compile the code: gcc -o testc99 main.c -O0

// FILE: main.c
// Description: Tests C99 IEEE 754 compliance.
// Note: use -O0 when compiling to make most of the tests pass
// Date: 11-24-2016

#include <stdio.h>
#include <math.h>
#include <inttypes.h>
#include <fenv.h>
#include <float.h>

#pragma STDC FENV_ACCESS ON

// Used to convert unsigned integer <--> double 
union Converter
{
    double d;
    uint64_t i;
};

typedef union Converter Converter;

#pragma STDC FENV_ACCESS on

int check_for_exceptions(void)
{
    // Set to 1 to enable debug printing
    if (0) {
        if (fetestexcept(FE_DIVBYZERO)) {
            printf("FE_DIVBYZERO detected\n");
        }
        if (fetestexcept(FE_INEXACT)) {
            printf("FE_INEXACT detected\n");
        }
        if (fetestexcept(FE_INVALID)) {
            printf("FE_INVALID detected\n");
        }
        if (fetestexcept(FE_OVERFLOW)) {
            printf("FE_OVERFLOW detected\n");
        }
        if (fetestexcept(FE_UNDERFLOW)) {
            printf("FE_UNDERFLOW detected\n");
        }
    }
    return fetestexcept(FE_ALL_EXCEPT);
}

// Adds two really big numbers together in order to cause an overflow
void test_overflow(void)
{
    double answer, num1, num2;
    num1 = 1.7 * pow(10, 308);  // the very limits of IEEE 754 double precision numbers
    num2 = num1;
    feclearexcept(FE_ALL_EXCEPT);
    // adding two really big numbers together should trigger an overflow
    answer = num1 + num2;
    printf("Test overflow...");
    if (check_for_exceptions() == (FE_OVERFLOW | FE_INEXACT)) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

void test_underflow(void)
{
    double answer;
    feclearexcept(FE_ALL_EXCEPT);
    //answer = DBL_MIN / 3000000000000000.0;    // does not produce an exception
    answer = fma(DBL_MIN, 1.0/10.0, 0);     // Inexact and underflow exceptions produced
    printf("Test underflow...");
    if (check_for_exceptions() == (FE_UNDERFLOW | FE_INEXACT)) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// Test if the inexact exception can be produced under the right conditions
void test_inexact(void)
{
    double answer;
    feclearexcept(FE_ALL_EXCEPT);
    answer = log(1.1);
    printf("Test inexact...");
    if (check_for_exceptions() == FE_INEXACT) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// test to see if the invalid exception works
void test_invalid(void)
{
    double d;
    feclearexcept(FE_ALL_EXCEPT);
    d = sqrt(-1.0);
    printf("Test invalid...");
    if (check_for_exceptions() == FE_INVALID) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// test the fused multiply-add operation
void test_fma(void)
{
    double result, correct_answer, num1, num2, num3;
    int iteration, max_iterations;

    result = 0.0;
    correct_answer = -13819435189200605973481192570224640.0;
    num1 = 5.2;
    num2 = 1.3 * pow(10, 18);
    num3 = -3.0;
    max_iterations = pow(10, 7);

    feclearexcept(FE_ALL_EXCEPT);

    // Test large number of multiplication, addition, and subtraction calculations
    for (iteration = 0; iteration < max_iterations; iteration++) {
        result += fma(num1, num2, num3);
        num1 += 1.0000002;
        num2 -= 1.3044 * pow(10,14);
        num3 += 0.953343;
    }

    // Test division - or multiplication with the reciprical
    result = fma(result, 1.0/3.14159265, -987654321.123456789);

    printf("Test fma()...");
    if (result == correct_answer) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// Test what happens with infinity - infinity
void test_inf_minus_inf()
{
    double answer;
    //answer = INFINITY - INFINITY;         // does not cause an exception, but does make a nan
    feclearexcept(FE_ALL_EXCEPT);
    answer = fma(INFINITY, 1, -INFINITY);   // does cause an invalid exception, answer is nan. -INFINITY - INFINITY doesn't cause an exception
    printf("Testing infinity - infinity...");
    if (check_for_exceptions() == FE_INVALID) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// Test signalling nan - should produce an invalid exception
void test_snan()
{
    double result;
    Converter c;
    c.i = 0x7FF0000000000001; 
    feclearexcept(FE_ALL_EXCEPT);
    result = fma(c.d, 10.4, 0.11);
    printf("Test snan...");
    if (check_for_exceptions() == FE_INVALID) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// Test quiet nan - no exceptions should be produced
void test_qnan()
{
    Converter c;
    double result;    
    c.i = 0x7fffffff;
    feclearexcept(FE_ALL_EXCEPT);
    result = fma(c.d, 10.4, 0.11);
    printf("Test qnan...");
    if (check_for_exceptions() == 0) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// Test infinity * zero for invalid exception
void test_inf_times_zero()
{
    double answer;
    //answer = INFINITY * 0;   // answer is nan, no exception
    feclearexcept(FE_ALL_EXCEPT);
    answer = fma(INFINITY, 0, 0); // answer is nan, invalid exception raised
    printf("Test infinity * 0...");
    if (check_for_exceptions() == FE_INVALID) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// Test division by zero exception
void test_one_divided_by_zero()
{
    double answer;
    feclearexcept(FE_ALL_EXCEPT);
    answer = fma(1.0, 1.0/0.0, 0.0);       // division by zero exception, answer = inf
    //answer = 1.0/0.0;                     // division by zero exception, answer = inf
    printf("Test division by zero...");
    if (check_for_exceptions() == FE_DIVBYZERO) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

// verify all rounding modes work correctly
void test_rounding_modes(void)
{
    double result, expected_result;

    printf("Test rounding...");
    do {
        fesetround(FE_TONEAREST);
        expected_result = 2.0;
        result = rint(2.1);        
        if (result != expected_result) {
            break;
        }

        fesetround(FE_UPWARD);
        expected_result = 3.0;
        result = rint(2.1);
        if (result != expected_result) {
            break;
        }

        fesetround(FE_DOWNWARD);
        expected_result = 2.0;
        result = rint(2.1);
        if (result != expected_result) {
            break;
        }

        fesetround(FE_TOWARDZERO);
        expected_result = 2.0;
        result = rint(2.1);
        if (result != expected_result) {
            break;
        }

        printf("pass\n");
        return;
    } while (0);
    printf("fail\n");
}

// Test square root results
void test_sqrt(void)
{
    double answer, result, base, expected_result;
    base = 8.123456;
    answer = pow(base, 2);
    result = sqrt(answer);
    expected_result = 8.1234559999999973;
    printf("Test sqrt...");
    if (result == expected_result) {
        printf("pass\n");
    } else {
        printf("fail\n");
    }
}

int main (int argc, const char * argv[]) {
    test_inf_minus_inf();
    test_inf_times_zero();
    test_one_divided_by_zero();
    test_fma();
    test_overflow();
    test_underflow();
    test_inexact();
    test_invalid();
    test_snan();
    test_qnan();
    test_rounding_modes();
    test_sqrt();
    return 0;
}

Answer 1

Note that Standard allows compilers to optimize floating-point computations across calls to fenv.h functions (see BZ 34678 for more details). To force compiler to respect them use -frounding-math and/or pragma STDC FENV ACCESS.

You can also try common tricks to hide constants from compiler (e.g. use volatile or clobber values via calls to dummy external function).

Answer 2

First, gcc does not support FENV_ACCESS (which you didn't write correctly btw, line 23 fails to compile because "on" has to use caps - why is it repeated anyway?)

Second, C allows fp expression results to vary with optimization levels because intermediate computations are both permitted to be fused and permitted to use larger size/precision (as in x87).

On Oracle, which is the most complete C99 compiler I have at hand - it even supports imaginary numbers! - I get fails on fma_test (32-bit x86 only) and qnan_test (32- and 64-bit sparc, 32- and 64-bit x86) and pass on the rest of your tests with every optimization level. The qnan test fails because it raises FE_INEXACT, the fma test on 32-bit x86 has large, but constant diff.

This is not a complete answer, but it was too large for a comment.. Perhaps an acceptable answer is: don't try, instead program in a way that is numerically stable and tolerant of the permitted variations.