c - 识别使 double 方程x + a == b成立的最小和最大x的最快算法

在静态分析的上下文中，我有兴趣在下面条件的then分支中确定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个浮点加法要便宜。

最佳答案

你已经在你的问题上给出了一个漂亮而优雅的解决方案：
如果用理性计算是便宜的，那就很简单：价值
因为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就不会在任何时候出现疯狂。
如果没有溢出和下溢，您可以从ulp1(b)获得b - nextafter(b, -1.0/0.0)。您可以从ulp2(b)获取nextafter(b, 1.0/0.0) - 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至少与b的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;
    }
  }
}