floating-point - 确定 C/C++ 中的浮点运算是否发生舍入

Question

我正在尝试提出一种有效的方法来确定 IEEE-754 操作何时会/确实会发生舍入。不幸的是我无法简单地检查硬件标志。它必须在几个不同的平台上运行。

我想到的方法之一是通过不同的舍入模式进行运算来比较结果。

添加示例:

    double result = operand1 + operand2;
    // save rounding mode
    int savedMode = fegetround();
    fesetround(FE_UPWARD);
    double upResult = operand1 + operand2;
    fesetround(FE_DOWNWARD);
    double downResult = operand1 + operand2;
    // restore rounding mode
    fesetround(savedMode);
    return (result != upResult) || (result != downResult);

但这显然效率很低，因为它必须执行 3 次操作。

Answer 1

您的示例不一定会通过优化给出正确的结果 级别 -O1 或更高。看这个Godbolt link : 编译器只生成一个附加vaddsd。

经过优化 level -O0 程序集看起来没问题，但这会导致代码效率低下。 而且调用fegetround和fesetround的成本相对较高， 与一些浮点运算的成本相比。

下面的( self 解释的)代码可能是一个有趣的选择。 它使用著名的算法2Sum和2ProdFMA。在没有硬件 fma 或 fma 仿真的系统上，您可以使用 2Prod 算法代替 2ProdFMA， 例如，请参阅精确浮点乘积和求幂， 作者:Stef Graillat。

/*
gcc -m64 -Wall -O3 -march=haswell round_ex.c -lm
   or with fma emulation on systems without hardware fma support, for example:
gcc -m64 -Wall -O3  -march=nehalem  round_ex.c -lm
*/

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

int add_is_not_exact(double operand1, double operand2){
    double a = operand1;
    double b = operand2;
    double s, t, a_1, b_1, d_a, d_b;
    /* Algorithm 2Sum computes s and t such that a + b = s + t, exactly.         */
    /* Here t is the error of the floating-point addition s = a + b.             */
    /* See, for example, On the robustness of the 2Sum and Fast2Sum algorithms   */
    /* by Boldo, Graillat, and Muller                                            */
    s = a + b;
    a_1 = s - b;
    b_1 = s - a_1;
    d_a = a - a_1;
    d_b = b - b_1;
    t = d_a + d_b;
    return (t!=0.0);
}


int sub_is_not_exact(double operand1, double operand2){
    return add_is_not_exact(operand1, -operand2);
}


int mul_is_not_exact(double operand1, double operand2){
    double a = operand1;
    double b = operand2;
    double s, t;
    /* Algorithm 2ProdFMA computes s and t such that a * b = s + t, exactly.     */
    /* Here t is the error of the floating-point multiplication s = a * b.       */
    /* See, for example, Accurate Floating Point Product and Exponentiation      */
    /* by Graillat                                                               */
    s = a * b;
    t = fma(a, b, -s);
    if (s!=0) return (t!=0.0);       /* No underflow of a*b                                */
    else return (a!=0.0)&&(b!=0.0);  /* Underflow: inexact if s=0, but (a!=0.0)&&(b!=0.0)  */
}


int div_is_not_exact(double operand1, double operand2){
    double a = operand1;
    double b = operand2;
    double s, t;
    s = a / b;
    t = fma(s, b, -a);  /* fma(x,y,z) computes x*y+z with infinite intermediate precision */
    return (t!=0.0);
}


int main(){

    printf("add_is_not_exact(10.0, 1.0) = %i\n", add_is_not_exact(10.0, 1.0));
    printf("sub_is_not_exact(10.0, 1.0) = %i\n", sub_is_not_exact(10.0, 1.0));
    printf("mul_is_not_exact( 2.5, 2.5) = %i\n", mul_is_not_exact( 2.5, 2.5));
    printf("div_is_not_exact(  10, 2.5) = %i\n", div_is_not_exact(  10, 2.5));
    printf("add_is_not_exact(10.0, 0.1) = %i\n", add_is_not_exact(10.0, 0.1));
    printf("sub_is_not_exact(10.0, 0.1) = %i\n", sub_is_not_exact(10.0, 0.1));
    printf("mul_is_not_exact( 2.6, 2.6) = %i\n", mul_is_not_exact( 2.6, 2.6));
    printf("div_is_not_exact(  10, 2.6) = %i\n", div_is_not_exact(  10, 2.6));

    printf("\n0x1.0p-300 = %20e, 0x1.0p-600 = %20e \n", 0x1.0p-300 , 0x1.0p-600 );
    printf("mul_is_not_exact( 0x1.0p-300, 0x1.0p-300) = %i\n", mul_is_not_exact( 0x1.0p-300, 0x1.0p-300));
    printf("mul_is_not_exact( 0x1.0p-600, 0x1.0p-600) = %i\n", mul_is_not_exact( 0x1.0p-600, 0x1.0p-600));

}

输出为:

$ ./a.out
add_is_not_exact(10.0, 1.0) = 0
sub_is_not_exact(10.0, 1.0) = 0
mul_is_not_exact( 2.5, 2.5) = 0
div_is_not_exact(  10, 2.5) = 0
add_is_not_exact(10.0, 0.1) = 1
sub_is_not_exact(10.0, 0.1) = 1
mul_is_not_exact( 2.6, 2.6) = 1
div_is_not_exact(  10, 2.6) = 1

0x1.0p-300 =         4.909093e-91, 0x1.0p-600 =        2.409920e-181 
mul_is_not_exact( 0x1.0p-300, 0x1.0p-300) = 0
mul_is_not_exact( 0x1.0p-600, 0x1.0p-600) = 1

---

正如评论中所指出的，也可以直接阅读 控制和状态寄存器:

#include <fenv.h>
#pragma STDC FENV_ACCESS ON

int add_is_not_exact_v2(double a, double b)
{    
    fexcept_t excepts;
    feclearexcept(FE_ALL_EXCEPT);
    double c = a+b;
    int tst = fetestexcept(FE_INEXACT);
    return (tst!=0);
}

但请注意，这可能不适用于编译器优化级别 -O1 或更高级别。 在这种情况下，addsd 双加指令有时会被完全优化掉， 导致错误的结果。 例如，对于 gcc 8.2 gcc -m64 -O1 -march=nehalem:

add_is_not_exact_v2:
        sub     rsp, 8
        mov     edi, 61
        call    feclearexcept
        mov     edi, 32
        call    fetestexcept
        test    eax, eax
        setne   al
        movzx   eax, al
        add     rsp, 8
        ret

优化级别-O0，有2个函数调用，并且相对 扩展指令来修改控制和状态寄存器，这不一定是最有效的解决方案。