除1外，IEEE 754浮点数机器epsilon的高效计算方法问询

问题背景 #

给定一个可被IEEE 754表示的正数a，求最大的可表示数e，使得按IEEE 754浮点算术规则计算时a + e == a成立。已知当a=1时，常被关联到std::numeric_limits<double>::epsilon()，但需明确：这个epsilon是1.0到下一个可表示浮点数的差值（即ULP(1.0)），而满足1 + e == 1的最大e其实是它的一半。当前可用朴素二分法实现，但需要比二分法更快的方案，且排除直接操作浮点数二进制位的“位魔法”。

最优解法：利用标准库函数直接计算 #

无需迭代二分，直接基于IEEE 754浮点数特性结合标准库函数即可O(1)时间完成计算，步骤如下：

• 计算a的ULP值：ULP（Unit in the Last Place）是相邻两个可表示浮点数的差值，用std::nextafter(a, std::numeric_limits<double>::infinity())可直接获取a的下一个更大的可表示浮点数，两者的差值就是ULP(a)。
• 得到最大e：满足a + e == a的最大e为ULP(a)的一半，即ulp_a / 2.0。

代码实现 #

#include <limits>
#include <cmath>

double max_e_that_leaves_a_unchanged(double a) {
    const double next_higher = std::nextafter(a, std::numeric_limits<double>::infinity());
    const double ulp = next_higher - a;
    return ulp / 2.0;
}

原理说明 #

IEEE 754默认采用向最近偶数舍入规则：当e小于等于ULP(a)/2时，a + e的结果会被舍入回a；一旦e超过这个值，a + e就会被舍入到a的下一个可表示浮点数，此时a + e != a。而ULP(a)/2本身是可被IEEE 754表示的（std::nextafter也能正确处理非规格化数场景）。

这种方法开销极低，std::nextafter是编译器优化过的底层函数，执行速度远快于需要多次迭代的二分法。

内容的提问来源于stack exchange，提问作者Eli4ph