检查浮点除法是否以 +-inf 结尾的可移植方法

Question

我有一个形式1.0f / x 为xa的浮点除法float。我如何事先检查是否x非常接近0.0f结果将是 +-inf / undefined？我不确定标准限制中的 epsilon 是否足够。

问候。

Answer 1

先决条件

\n

C++ 不强制要求 IEEE-754 或特定的舍入方法。对于这个答案，我假设 IEEE-754 使用二进制格式并舍入到最接近的值，连到偶数。

\n

结论

\n

1/x当且仅当溢出fabs(x) <= std::ldexp(1, -std::numeric_limits<float>::max_exponent)。对于常量表达式，您可以使用std::numeric_limits<float>::min()/4.

\n

讨论

\n

在有限范围末尾选择舍入方向就好像指数范围继续下去一样。例如，使用十进制进行说明，如果最高可表示的有限数是 9.99\xe2\x80\xa210 ¹⁷，则下一个可表示的数（如果指数不受限制）将为 1.00\xe2\x80\xa210 ¹⁸。这两者之间的中点是 9.995\xe2\x80\xa210 ¹⁷，因此低于该值的数字向下舍入，高于该值的数字向上舍入。如果是偶数，则 9.995\xe2\x80\xa210 ¹⁷会向上舍入。

\n

对于二进制格式，最大可表示值为 (2\xe2\x88\x92\xce\xb5)\xe2\x80\xa22 ^q，其中 \xce\xb5 是 \xe2\x80\x9cmachine epsilon\xe2\x80\ x9d（1 的 ULP，因此 2-\xce\xb5 是最大可表示有效数），q是最大指数。那么发生舍入的点是 (2\xe2\x88\x92\xc2\xbd\xce\xb5)\xe2\x80\xa22 ^q。

\n

对于正x，如果 1/ x < (2\xe2\x88\x92\xc2\xbd\xce\xb5)\xe2\x80\xa22 ^q，结果向下舍入。否则，向上舍入为 \xe2\x88\x9e。因此，结果小于 \xe2\x88\x9e iff x > 1/((2\xe2\x88\x92\xc2\xbd\xce\xb5)\xe2\x80\xa22 ^q ) = 2 ^{\xe2\x88 \x92 q} /(2-\xc2\xbd\xce\xb5)。

\n

1/(2-\xc2\xbd\xce\xb5) 略大于\xc2\xbd，小于\xc2\xbd\xce\xb5，因此小于或等于它的最大可表示值为\xc2\ xbd。因此，对于正x， 的结果1/x小于 \xe2\x88\x9e iff x > 2 ^{\xe2\x88\x92 q} /2 = 2 ^{\xe2\x88\x92 q \xe2\x88\x921}，并且情况对于负x是对称的。

\n

C++ 告诉我们最大指数std::numeric_limits<double>::max_exponent（在标头中定义<limits>）。但是，C++ 将此指数校准为有效数范围 [\xc2\xbd, 1) 而不是 IEEE-754\xe2\x80\x99s [1, 2)，因此它比 q 大1。因此我们想要的 \xe2\x88\x92 q \xe2\x88\x921 就是-std::numeric_limits<double>::max_exponent。

\n

我们可以使用函数（在 中声明）计算 2 ^{\xe2\x88\x92 q \xe2\x88\x921}。ldexp<cmath>std::ldexp(1, -std::numeric_limits<float>::max_exponent)

\n

在 Apple Clang 11 中，该程序：

\n

#include <cmath>\n#include <iomanip>\n#include <iostream>\n#include <limits>\n\n\nint main(void)\n{\n    float x = std::ldexp(1, -std::numeric_limits<float>::max_exponent);\n\n    std::cout << std::setprecision(20) << x << " is too small, result will overflow:\\n";\n    std::cout << "\\t" << 1/x << ".\\n";\n\n    x = std::nexttoward(x, INFINITY);\n\n    std::cout << std::setprecision(20) << x << " is just big enough, result will not overflow:\\n";\n    std::cout << "\\t" << 1/x << ".\\n";\n}\n

\n

产生：

\n

\n2.9387358770557187699e-39 太小，结果会溢出:\n inf。\n2.9387372783541830947e-39 足够大，结果不会溢出:\n 3.4028220466166163425e+38。\n

\n

考虑到负数也会1/x溢出 iff fabs(x) <= std::ldexp(1, -std::numeric_limits<float>::max_exponent)。

\n

由于 IEEE-754 指定指数范围的方式，std::ldexp(1, -std::numeric_limits<float>::max_exponent)等于std::numeric_limits<float>::min()/4。（IEEE-754 指定最小正规指数为 1\xe2\x88\x92 q，因此我们希望的 \xe2\x88\x92 q \xe2\x88\x921 为 (1- q )-2。）

\n

Answer 2

我们可以通过反复试验来搜索极限：

#include <iostream>
#include <limits>

#include <cmath>

int main() {
    float limit = 0.0f;
    float result = 1.0f / limit;
    while (
        result == std::numeric_limits<float>::infinity()
        or std::isnan(result)
    ) {
        limit = std::nextafter(limit, 1.0f);
        result = 1.0f / limit;
    }
    std::cout << "Limit = " << limit << std::endl;
    std::cout << "1.0f / Limit = " << 1.0f / limit << std::endl;
}

这在我的系统上输出：

Limit = 2.93874e-39
1.0f / Limit = 3.40282e+38

然而，这并不是一个非常有效的解决方案。如果我们能够制定这个算法constexpr，这将缓解这个问题，但不幸的std::nextafter()是事实并非如此constexpr。

如果您知道您的环境正在使用 IEEE-754 空气算法，这些限制可能是恒定的，但当您要求可移植性时，我们不能总是这样假设。