c - C99 浮動小数点コードが最適化レベルで結果を変更しないようにする方法

Question

システムが C99 の IEEE 754 に準拠しているかどうかをテストするために、以下のコードを作成しました。問題は、最適化レベルが異なると結果が変わることです。最適化レベル 0 に固執すると、ほとんどのテストに合格します。しかし、最適化レベルを -O3 に上げると、ほとんどのテストが失敗します。-float-store を試しました。結果は変わりませんでした。

GCC 5.2.0 を使用する Mac OS 10.4.11 PowerPC と、GCC 4.2.1 および GCC 4.9.2 を使用する Mac OS 10.6.8 でこのプログラムを試しました。結果は同じです。ほとんどのテストは、最適化レベル 0 (-O0) を使用した場合にのみパスします。これは GCC のバグだと思いますが、報告する前に、コードに関する他の意見を聞きたいと思います。

これは私がコードをコンパイルする方法です: 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;
}