1

我在下面编写了这段代码来测试系统是否符合 C99 的 IEEE 754。问题是结果会随着优化级别的不同而变化。如果我坚持优化级别为零,则大多数测试都会通过。但是如果我将优化级别提高到 -O3,那么大多数测试都会失败。我试过-ffloat-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;
}
4

2 回答 2

1

首先,gcc 不支持 FENV_ACCESS (顺便说一句,您没有正确编写它,第 23 行无法编译,因为“on”必须使用大写字母 - 为什么还是要重复?)

其次,C 允许 fp 表达式结果随优化级别而变化,因为中间计算既允许融合也允许使用更大的尺寸/精度(如 x87 中)。

在 Oracle 上,这是我手头最完整的 C99 编译器——它甚至支持虚数!- 我在 fma_test(仅限 32 位 x86)和 qnan_test(32 位和 64 位 sparc,32 位和 64 位 x86)上失败,并通过每个优化级别的其余测试。qnan 测试失败,因为它引发了 FE_INEXACT,32 位 x86 上的 fma 测试具有很大但恒定的差异。

这不是一个完整的答案,但它对于评论来说太大了。也许可以接受的答案是:不要尝试,而是以数字稳定和容忍允许的变化的方式编程。

于 2016-12-30T17:23:28.357 回答
1

请注意,Standard 允许编译器在调用 fenv.h 函数时优化浮点计算(有关详细信息,请参阅BZ 34678 )。强制编译器尊重他们使用-frounding-math和/或编译指示 STDC FENV ACCESS

您还可以尝试常用技巧来隐藏编译器中的常量(例如,通过调用虚拟外部函数来使用 volatile 或 clobber 值)。

于 2016-12-26T19:57:29.017 回答