如何确定是在四舍五入前还是四舍五入后检测到微小或无法确定？

问题描述

IEEE 754-2008：

7.5 下溢

当检测到微小的非零结果时，应发出下溢异常信号。对于二进制格式，这应该是：

a) 舍入后 — 当计算出非零结果时，就好像指数范围是无界的一样 将严格介于 ±b^emin 或

b) 舍入前 — 当计算出非零结果时，就好像指数范围和 无界精度将严格介于 ±b^emin 之间。

实现者应选择检测微小的方式，但应以相同的方式检测基数为 2 的所有操作的微小，包括二进制舍入属性下的转换操作。

但是，C11 和 C17..C2x（工作草案 - 2020 年 2 月 5 日，n2479.pdf）都没有提及微小：

$ pdfgrep.exe -i 'tininess' ISO-IEC-9899-2011.pdf n2479.pdf --color never
<nothing>

困惑。

问题：为什么没有 FLT_tininesS 宏 (-1 -- indeterminable,0 -- after rounding,1 -- before rounding)？

更新。问题原因：像往常一样：一些 FP 测试（测试 FP 操作生成的结果的正确性）失败，因为预期引发的异常是 FE_INEXACT，而实际引发的异常是 FE_INEXACT 和 {{1 }}。然后结果是 HW FE_UNDERFLOW。因此，逻辑问题出现了：“如何确定是在四舍五入之前还是在四舍五入之后或不可确定的？”。由于在编译时无法确定，所以需要在运行时确定。

解决方法

以下程序可以确定是在四舍五入之前还是之后报告微小。

#include <fenv.h>
#include <float.h>


_Static_assert(FLT_RADIX == 2,"This program expects binary floating-point.");

typedef float Float;
enum {  //  Change FLT prefix according to type set for Float,above.
    Precision       = FLT_MANT_DIG,// Number of bits in significand.
    MinimumExponent = FLT_MIN_EXP-1,// Minimum normal exponent.
        /*  The -1 is due to C's definition of floating-point exponents being
            for significands in [1/2,1) instead of [1,2).
        */
};


// Use the following if your compiler supports it.  Not all do.
//#pragma STDC FENV_ACCESS ON


//  Report true iff a*b reports underflow.
static _Bool ProductUnderflows(Float a,Float b)
{
    feclearexcept(FE_ALL_EXCEPT);
    volatile Float c;
    c = a*b;
    return fetestexcept(FE_UNDERFLOW);
}


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


int main(void)
{
    if (fesetround(FE_TONEAREST) != 0)
    {
        fprintf(stderr,"Error,cannot set rounding mode to nearest.\n");
        exit(EXIT_FAILURE);
    }

    /*  Find the least positive integer that does not divide the number of bits
        in a significand (also called p or the precision of the type).
    */
    int q = 1;
    while (Precision % q == 0)
        ++q;

    //  Set a to a string of q bits after the radix point.
    Float a = 1 - ldexp(1,-q);

    /*  Consider 1/a.  This necessarily rounds down and sets b to a repeating
        pattern of a 1 bit followed by q-1 0 bits.

        To see that it rounds down,consider the binary representation of the
        mathematical quotient 1/a.  It is a repeating pattern of a 1 bit
        followed by q-1 0 bits.  So the 1 bits land at offsets from the first 1
        bit of q,2q,3q,and so on.  So they only land at multiples of q.  And
        we know p is not a multiple of q,so there is no 1 bit at the position
        p bits beyond the leading bit.  In other words,the first bit that is
        does not fit in the p-bit significand is 0. So the residue being
        discarded during rounding is less than 1/2 ULP,so round-to-nearest
        rounds down.

        We set b to 1/a scaled so that a*b is just below the normal range.

        Then the mathematical product of a and b has a significand of
        ceil(p/q)*q 1 bits,which is greater than p,so the product must be
        rounded to fit in a signifcand.  In round-to-nearest-ties-to-even mode,it will round upward,so the floating-point product of a and b will be
        the smallest normal number.  Therefore,there is an underflow if
        tininess is detected before rounding but not if it is detected after
        rounding.
    */
    Float b = ldexp(1/a,MinimumExponent);

    printf("a = %a.\n",a);
    printf("b = %a.\n",b);

    /*  Test that we hit the boundary correctly:  (a/2)*b underflows but
        (2*a)*b does not.  Also test that underflow reporting works.
    */
    if (!ProductUnderflows(a/2,b))
    {
        fprintf(stderr,"Internal error,%a * %a -> %a is expected to underflow but did not.\n",a/2,b,(a/2)*b);
        exit(EXIT_FAILURE);
    }
    if (ProductUnderflows(2*a,%a * %a -> %a is expected not to underflow but did.\n",2*a,(2*a)*b);
        exit(EXIT_FAILURE);
    }

    //  Test whether tininess is detected before or after rounding.
    printf("Tininess is detected %s rounding.\n",ProductUnderflows(a,b) ? "before" : "after");
}