检查浮点除法是否以+-inf结尾的便携式方法

Question

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

问候。

Answer 1

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

#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 airthmetic，那么这些限制可能是不变的，但是当您要求可移植性时，我们不能总是这样假设。

Answer 2

先决条件

C++ 不强制要求 IEEE-754 或特定的舍入方法。对于这个答案，我假设 IEEE-754 与二进制格式和四舍五入到最近的平局一起使用。

结论

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

讨论

在有限范围的末尾进行舍入，就好像指数继续前进一样。例如，使用十进制来说明，如果最高可表示的有限数是 9.99•10¹⁷，那么如果指数不受限制，则下一个可表示的数是 1.00•10¹⁸。这两者之间的中点是 9.995•10¹⁷，因此低于该数字的数字向下舍入，高于该数字的数字向上舍入。平局时，9.995•10¹⁷ 向上取整。

对于二进制格式，最大可表示值为 (2−ε)•2^q，其中 ε 是“机器 epsilon”（1 的 ULP，所以 2-ε 是最大的可表示有效数）并且 q 是最大指数。那么发生舍入的点是 (2−½ε)•2^q。

如果 1/x q，则结果向下舍入。否则，向上舍入到∞。因此，结果小于 ∞ 当且仅当 x > 1/((2−½ε)•2^q) = 2^{− q}/(2-½ε).

1/(2-½ε) 略大于½，小于½ε，因此小于或等于它的最大可表示值为½。因此，1/x 的结果小于 ∞ 当且仅当 x > 2^-q/2 = 2^{- q-1}.

C++ 用std::numeric_limits<double>::max_exponent 告诉我们最大指数（在标题<limits> 中定义）。然而，C++ 将这个指数校准为 [½, 1) 的有效数字范围，而不是 IEEE-754 的 [1, 2)，因此它比 q 大一。因此我们想要的-q-1 就是-std::numeric_limits<double>::max_exponent。

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

使用 Apple Clang 11，此程序：

#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>


int main(void)
{
    float x = std::ldexp(1, -std::numeric_limits<float>::max_exponent);

    std::cout << std::setprecision(20) << x << " is too small, result will overflow:\n";
    std::cout << "\t" << 1/x << ".\n";

    x = std::nexttoward(x, INFINITY);

    std::cout << std::setprecision(20) << x << " is just big enough, result will not overflow:\n";
    std::cout << "\t" << 1/x << ".\n";
}

产生：

2.9387358770557187699e-39 太小，结果会溢出： 信息。 2.9387372783541830947e-39 刚好够大，结果不会溢出： 3.4028220466166163425e+38。

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

由于 IEEE-754 指定指数范围的方式，std::ldexp(1, -std::numeric_limits<float>::max_exponent) 等于 std::numeric_limits<float>::min()/4。 （IEEE-754 规定最小正态指数为 1-q，所以我们想要的 -q-1 是 (1-q) -2.)

Answer 3

由于您正在寻找常数值，我们实际上可以使用 SMT 求解器来找到 x 的最小/最大值，其中除法 1/x 将产生无穷大。我使用 Microsoft 的 z3 SMT 求解器 (https://github.com/Z3Prover/z3) 完成了这项工作，并从 Haskell SBV 绑定 (http://leventerkok.github.io/sbv/) 编写了脚本。我明白了：

Prelude Data.SBV> optimize Lexicographic $ do x <- sFloat "x"; constrain (fpIsInfinite (1/x)); minimize "x" x
Optimal model:
  x   = -2.938736e-39 :: Float
  x_0 =    2145386495 :: Word32
Prelude Data.SBV> optimize Lexicographic $ do x <- sFloat "x"; constrain (fpIsInfinite (1/x)); maximize "x" x
Optimal model:
  x   = 2.938736e-39 :: Float
  x_0 =   2149580800 :: Word32

如果你眯着眼睛看，你会发现它建议的值介于-2.938736e-39 和2.938736e-39 之间，其中1/x 变为无穷大。

如果你想写这个“可移植”而没有任何舍入问题，你应该使用十六进制表示法，即-0x1p-128和0x1p-128。

我相信这些数字符合@Eric Postpischill 的价值观；他的分析当然非常有用，但这是使用自动定理证明技术找到此类值的另一种方法。