Как предотвратить код с плавающей запятой C99 от изменения результатов с уровнем оптимизации

Question

Я сделал этот код ниже, чтобы проверить соответствие системы CEE IEEE 754. Проблема заключается в изменении результатов с разными уровнями оптимизации. Если я придерживаюсь нулевого уровня оптимизации, большинство тестов проходят. Но если я повышу уровень оптимизации до -O3, то большинство тестов не удастся. Я попробовал -ффлоат-магазин. Это не изменило никаких результатов.

Я пробовал эту программу на 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.

Вы также можете попробовать общие трюки, чтобы скрыть константы из компилятора (например, использовать изменчивые или clobber-значения посредством вызовов фиктивной внешней функции).

Answer 2

Во-первых, НКУ не поддерживает FENV_ACCESS (который вы не писали правильно кстати, строка 23 не может компилироваться, так как «на» должен использовать заглавные буквы - почему повторяется она так или иначе)

Во-вторых, C позволяет Результаты выражения fp варьируются в зависимости от уровней оптимизации, потому что промежуточным вычислениям разрешено сливаться и разрешено использовать больший размер/точность (как в x87).

В Oracle, который является самым полным компилятором C99, который у меня под рукой, он поддерживает даже мнимые числа! - Я терпит неудачу в fma_test (только 32-разрядный x86) и qnan_test (32- и 64-разрядный sparc, 32- и 64-разрядный x86) и передаю остальные ваши тесты с каждым уровнем оптимизации. Тест qnan терпит неудачу, поскольку он вызывает FE_INEXACT, тест fma на 32-разрядном x86 имеет большой, но постоянный diff.

Это не полный ответ, но он слишком велик для комментария. Возможно, приемлемый ответ: не пытайтесь, а программируйте таким образом, чтобы он был численно устойчивым и терпимым к допустимым вариациям.