c - 各种浮点值的 epsilon

Question

有一个最接近零的FLT_MIN常量。如何最接近某个数字值？

举个例子:

float nearest_to_1000 = 1000.0f + epsilon;
// epsilon must be the smallest value satisfying condition:
// nearest_to_1000 > 1000.0f

我更喜欢不使用特殊函数的数字公式。

Answer 1

C为此提供了一个函数，在<math.h>中 header 。 nextafterf(x, INFINITY)是 x 之后的下一个可表示值，朝 INFINITY 的方向.

但是，如果您愿意自己动手​​:

以下内容返回您寻求的 epsilon，对于单精度(浮点)，假设 IEEE 754。请参阅底部有关使用库例程的注释。

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


/*  Return the ULP of q.

    This was inspired by Algorithm 3.5 in Siegfried M. Rump, Takeshi Ogita, and
    Shin'ichi Oishi, "Accurate Floating-Point Summation", _Technical Report
    05.12_, Faculty for Information and Communication Sciences, Hamburg
    University of Technology, November 13, 2005.
*/
float ULP(float q)
{
    // SmallestPositive is the smallest positive floating-point number.
    static const float SmallestPositive = FLT_EPSILON * FLT_MIN;

    /*  Scale is .75 ULP, so multiplying it by any significand in [1, 2) yields
        something in [.75 ULP, 1.5 ULP) (even with rounding).
    */
    static const float Scale = 0.75f * FLT_EPSILON;

    q = fabsf(q);

    /*  In fmaf(q, -Scale, q), we subtract q*Scale from q, and q*Scale is
        something more than .5 ULP but less than 1.5 ULP.  That must produce q
        - 1 ULP.  Then we subtract that from q, so we get 1 ULP.

        The significand 1 is of particular interest.  We subtract .75 ULP from
        q, which is midway between the greatest two floating-point numbers less
        than q.  Since we round to even, the lesser one is selected, which is
        less than q by 1 ULP of q, although 2 ULP of itself.
    */
    return fmaxf(SmallestPositive, q - fmaf(q, -Scale, q));
}

下面的代码返回传递的值之后可以用 float 表示的下一个值(将 -0 和 +0 视为相同)。

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


/*  Return the next floating-point value after the finite value q.

    This was inspired by Algorithm 3.5 in Siegfried M. Rump, Takeshi Ogita, and
    Shin'ichi Oishi, "Accurate Floating-Point Summation", _Technical Report
    05.12_, Faculty for Information and Communication Sciences, Hamburg
    University of Technology, November 13, 2005.
*/
float NextAfterf(float q)
{
    /*  Scale is .625 ULP, so multiplying it by any significand in [1, 2)
        yields something in [.625 ULP, 1.25 ULP].
    */
    static const float Scale = 0.625f * FLT_EPSILON;

    /*  Either of the following may be used, according to preference and
        performance characteristics.  In either case, use a fused multiply-add
        (fmaf) to add to q a number that is in [.625 ULP, 1.25 ULP].  When this
        is rounded to the floating-point format, it must produce the next
        number after q.
    */
#if 0
    // SmallestPositive is the smallest positive floating-point number.
    static const float SmallestPositive = FLT_EPSILON * FLT_MIN;

    if (fabsf(q) < 2*FLT_MIN)
        return q + SmallestPositive;

    return fmaf(fabsf(q), Scale, q);
#else
    return fmaf(fmaxf(fabsf(q), FLT_MIN), Scale, q);
#endif
}

使用了库例程，但是 fmaxf (其参数的最大值)和 fabsf (绝对值)很容易被替换。 fmaf应该在具有融合乘加的架构上编译为硬件指令。如果失败的话，fmaf(a, b, c)在此使用中可以替换为 (double) a * b + c 。 (IEEE-754 二进制 64 具有足够的范围和精度来替换 fmaf 。 double 的其他选择可能没有。)

融合乘加的另一种替代方法是为 q * Scale 的情况添加一些测试。将是次正常的并单独处理它们。其他情况可以用普通*单独进行乘法和加法。和+运算符。