识别使双精度等式成立的最小和最大x的最快算法x + a == b true

Question

在静态分析的环境中，我有兴趣确定条件的当时分支中的x的值下面：

double x; 
x = …; 
if (x + a == b) 
{ 
    … 

a和b可以假定为双精度常数（一般化到任意的表达式是问题的最容易的部分），并且可以假定编译器遵循IEEE 754严格（FLT_EVAL_METHOD是0 ）。运行时的舍入模式可以假定为最接近偶数。如果使用有理数进行计算很便宜，那将很简单：x的值将是包含在有理区间中的双精度数（b - a - 0.5 * ulp1（b）... b - a + 0.5 * ulp2（b））。如果b是偶数，则应该包含边界，如果b是奇数，则排除边界，并且ulp1和ulp2是“ULP”的两个略有不同的定义，如果不介意在两个幂次方面失去一点精度，则可以采用相同的定义。

不幸的是，有理数计算可能会很昂贵。考虑另一种可能性是通过64位双精度加法（每个操作决定结果的一位）通过二分法获得每个边界。获得下限和上限的128个浮点加法可能比基于数学的任何解决方案都快。

我想知道是否有一种方法来改善“128点浮点加法”的想法。其实我有我自己的解决方案，涉及到四舍五入模式和nextafter电话的变化，但我不想抽象任何人的风格，并导致他们错过比我目前更优雅的解决方案。另外我不确定两次更改舍入模式实际上比64个浮点加法更便宜。

Answer 1

你已经给你的问题一个很好的和优雅的解决方案：

如果计算有理数很便宜，这将是简单的：值 对于x将包含在合理的双精度数（b-a-0.5 * ulp1（b）... b-a + 0.5 * ulp2（b））。如果b是偶数，则应该包括边界 ，如果b是奇数，则应该包括边界 ，并且ulp1和 ulp2是两个稍微不同的“ULP”的定义，如果不介意 的功率失去一点精度二。

以下是基于本段落的问题的部分解决方案的半理由草图。希望我有机会尽快充实它。要得到真正的解决方案，你必须处理低于正常值，零点，NaN和所有其他有趣的东西。我假设a和b是，例如，1e-300 < |a| < 1e300和1e-300 < |b| < 1e300，以便在任何点都不会发生疯狂。

没有上溢和下溢，您可以从b - nextafter(b, -1.0/0.0)得到ulp1(b)。你可以从nextafter(b, 1.0/0.0) - b得到ulp2(b)。

如果b/2 <= a <= 2b，那么Sterbenz定理告诉你b - a是确切的。所以(b - a) - ulp1/2将是最接近的double到下限和(b - a) + ulp2/2将是最接近的double到上限。尝试使用这些值，以及之前和之后的值，并选择最宽的工作时间间隔。

如果b > 2a,b - a > b/2。计算出的值b - a最多不超过半个ulp。一个ulp1是最多两个ulp，因为是一个ulp2，所以您给出的合理间隔最多为两个ulp。找出与b-a工作的五个最接近的值中的哪一个。

如果是a > 2b，那么b-a的ulp至少与ulpb的ulp一样大;如果有的话，我敢打赌它必须是b-a的三个最接近的值之一。我想象a和b具有不同的标志作品的情况相似。

我写了一小堆实现这个想法的C++代码。在无聊等待之前，它没有失败随机模糊测试（在几个不同的范围内）。这里是：

void addeq_range(double a, double b, double &xlo, double &xhi) { 
    if (a != a) return; // empty interval 
    if (b != b) { 
    if (a-a != 0) { xlo = xhi = -a; return; } 
    else return; // empty interval 
    } 
    if (b-b != 0) { 
    // TODO: handle me. 
    } 

    // b is now guaranteed to be finite. 
    if (a-a != 0) return; // empty interval 

    if (b < 0) { 
    addeq_range(-a, -b, xlo, xhi); 
    xlo = -xlo; 
    xhi = -xhi; 
    return; 
    } 

    // b is now guaranteed to be zero or positive finite and a is finite. 
    if (a >= b/2 && a <= 2*b) { 
    double upulp = nextafter(b, 1.0/0.0) - b; 
    double downulp = b - nextafter(b, -1.0/0.0); 
    xlo = (b-a) - downulp/2; 
    xhi = (b-a) + upulp/2; 
    if (xlo + a == b) { 
     xlo = nextafter(xlo, -1.0/0.0); 
     if (xlo + a != b) xlo = nextafter(xlo, 1.0/0.0); 
    } else xlo = nextafter(xlo, 1.0/0.0); 
    if (xhi + a == b) { 
     xhi = nextafter(xhi, 1.0/0.0); 
     if (xhi + a != b) xhi = nextafter(xhi, -1.0/0.0); 
    } else xhi = nextafter(xhi, -1.0/0.0); 
    } else { 
    double xmid = b-a; 
    if (xmid + a < b) { 
     xhi = xlo = nextafter(xmid, 1.0/0.0); 
     if (xhi + a != b) xhi = xmid; 
    } else if (xmid + a == b) { 
     xlo = nextafter(xmid, -1.0/0.0); 
     xhi = nextafter(xmid, 1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
     if (xhi + a != b) xhi = xmid; 
    } else { 
     xlo = xhi = nextafter(xmid, -1.0/0.0); 
     if (xlo + a != b) xlo = xmid; 
    } 
    } 
}