3

我有一个浮点除法，格式1.0f / x 为. 我将如何事先检查是否如此接近结果将是 +-inf / undefined？我不确定标准限制中的 epsilon 是否足够。xfloatx0.0f

问候。

4

3 回答 3

6

先决条件

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 17，那么如果指数不受限制，则下一个可表示的数将是1.00•10 18。这两者之间的中点是 9.995•10 17，因此低于的数字向下舍入，高于的数字向上舍入。平局时，9.995•10 17向上取整。

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

如果 1/ x < (2−½ε)•2 q，则结果向下舍入。否则，向上舍入到∞。因此，如果x > 1/((2−½ε)•2 q ) = 2 −<em>q /(2-½ε) ，则结果小于 ∞ 。

1/(2-½ε) 略大于½，小于½ε，因此小于或等于它的最大可表示值是½。因此，如果x > 2 -<em>q /2 = 2 -<em>q-11/x的结果小于 ∞ 。

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

我们可以用函数（在 中声明）计算 2 −<em>q−1：。ldexp<cmath>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−<em>q，因此我们想要的 −<em>q−1 为 (1- q )-2。）

于 2021-09-24T12:41:08.013 回答
6

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

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

于 2021-09-24T11:14:13.943 回答
0

x由于您正在寻找常数值，我们实际上可以使用 SMT 求解器来找到除法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-391/x

如果您想在没有任何舍入问题的情况下“可移植”地编写此代码，则应使用十六进制表示法，即-0x1p-128和0x1p-128.

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

于 2021-09-24T18:37:05.447 回答