如何防止C99浮點代碼通過優化級別改變結果

Question

我在下面做了下面的代碼來測試系統是否符合C99的IEEE 754.問題是結果隨着不同的優化級別而改變。如果我堅持零級優化，大部分測試都會通過。但如果我將優化級別提高到-O3，那麼大多數測試都會失敗。我試過了-ffloat-store。它沒有改變任何結果。

我已經在Mac OS 10.4.11 PowerPC上使用GCC 5.2.0和Mac OS 10.6.8使用GCC 4.2.1和GCC 4.9.2試過這個程序。結果是一樣的。大多數測試僅在使用優化級別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; 
} 

Answer 1

注意，標準允許編譯器在調用優化浮點計算到fenv.h功能（參見BZ 34678爲更多細節）。要強制編譯器尊重它們，請使用-frounding-math和/或pragma STDC FENV ACCESS。

您也可以嘗試使用常用技巧來隱藏編譯器中的常量（例如，通過調用虛擬外部函數來使用volatile或clobber值）。

Answer 2

首先，GCC不支持FENV_ACCESS（你沒有寫正確順便說一句，23行編譯失敗，因爲「上」必須使用帽 - ？爲什麼會重複反正）

其次，C允許fp表達式的結果隨優化級別而變化，因爲中間計算都允許被融合並允許使用更大的大小/精度（如x87）。

在Oracle上，這是我手頭最完整的C99編譯器 - 它甚至支持虛數！ - 我在fma_test（僅限32位x86）和qnan_test（32位和64位sparc，32位和64位x86）上失敗，並將每個優化級別的測試傳遞給其餘部分。 qnan測試失敗是因爲它引發FE_INEXACT，32位x86上的fma測試具有較大但恆定的差異。

這不是一個完整的答案，但它太大了評論..也許一個可以接受的答案是：不要嘗試，而是以一種在數值上穩定並容許允許變化的方式進行編程。