最快的算法来识别的最小和最大X中使双precision方程x + A == B真 [英] Fastest algorithm to identify the smallest and largest x that make the double-precision equation x + a == b true

问题描述

在静态分析的背景下，我感兴趣的是确定的值x 在以下条件的再分支：

In the context of static analysis, I am interested in determining the values of x in the then-branch of the conditional below:

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

A 和 B 可以被认为是双precision常数（推广到任意的前pressions是问题的最容易的部分），并且可以假设编译器遵循IEEE 754严格（ FLT_EVAL_METHOD 0）。在运行时的舍入模式可被假定为至最近，甚至

a and b can be assumed to be double-precision constants (generalizing to arbitrary expressions is the easiest part of the problem), and the compiler can be assumed to follow IEEE 754 strictly (FLT_EVAL_METHOD is 0). The rounding mode at run-time can be assumed to be to nearest-even.

如果计算有理数很便宜，这将是简单： X 的值将包含在合理的区间（双precision数b - 一个 - 0.5 * ULP1（二）。b - A + 0.5 * ulp2（b））的。该范围应包括如 B 均匀，排除在外，如果 B 是奇数，而且ULP1和ulp2是两个略有不同ULP可以采取相同的定义，如果一个人不介意失去对两个大国一点点precision。

If computing with rationals was cheap, it would be simple: the values for x would be the double-precision numbers contained in the rational interval (b - a - 0.5 * ulp1(b) … b - a + 0.5 * ulp2(b)). The bounds should be included if b is even, excluded if b is odd, and ulp1 and ulp2 are two slightly different definitions of "ULP" that can be taken identical if one does not mind losing a little precision on powers of two.

不幸的是，与有理数计算可能是昂贵的。考虑到的另一种可能性是通过二分法来获得每个界，在64双precision添加剂（每个操作判定的结果中的一个位）。 128浮点加法以获得下限和上限很可能比基于数学的任何解决方案更快。

Unfortunately, computing with rationals can be expensive. Consider that another possibility is to obtain each of the bounds by dichotomy, in 64 double-precision additions (each operation deciding one bit of the result). 128 floating-point additions to obtain the lower and upper bounds may well be faster than any solution based on maths.

我想知道是否有改善在128浮点加法的想法的一种方式。其实我有我的，涉及舍入模式和函数nextafter 呼吁改变自己的解决方案，但我不想抽筋任何人的风格，导致他们错过一个比一个更优雅的解决方案我现在有。另外，我不知道，改变舍入模式的两倍实际上价格比64浮点加法。

I am wondering if there is a way to improve over the "128 floating-point additions" idea. Actually I have my own solution involving changes of rounding mode and nextafter calls, but I wouldn't want to cramp anyone's style and cause them to miss a more elegant solution than the one I currently have. Also I am not sure that changing the rounding mode twice is actually cheaper than 64 floating-point additions.

推荐答案

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

You already gave a nice and elegant solution in your question:

如果计算有理数很便宜，这将是简单的值
  对于x将是包含在合理的双链precision号码
  间隔（二 - 一 - 0.5 * ULP1（二）。B - A + 0.5 * ulp2（b））的。该范围
  应包括如果B是偶数，排除在外，如果b为奇数，且ULP1和
  ulp2两种可以采取ULP的稍微不同的定义
  同样，如果不介意的权力失去了一点点precision
  两项。

If computing with rationals was cheap, it would be simple: the values for x would be the double-precision numbers contained in the rational interval (b - a - 0.5 * ulp1(b) … b - a + 0.5 * ulp2(b)). The bounds should be included if b is even, excluded if b is odd, and ulp1 and ulp2 are two slightly different definitions of "ULP" that can be taken identical if one does not mind losing a little precision on powers of two.

以下是部分解决了基于本款问题的半理性的草图。希望我得到一个机会，很快就割肉出来。要获得真正的解决方案，你必须处理次归，零，NaN的，和所有其他有趣的东西。我要去假设 A 和 B 是，比方说，这样 1e- 300℃ | A | ＆LT; 1e300 和 1E-300℃ | B | ＆LT; 1e300 所以没有疯狂出现在任何点。

What follows is a half-reasoned sketch of a partial solution to the problem based on this paragraph. Hopefully I'll get a chance to flesh it out soon. To get a real solution, you'll have to handle subnormals, zeroes, NaNs, and all that other fun stuff. I'm going to assume that a and b are, say, such that 1e-300 < |a| < 1e300 and 1e-300 < |b| < 1e300 so that no craziness occurs at any point.

缺席溢和下溢，即可获得 ULP1（B）从 B - 函数nextafter（B，-1.0 / 0.0）。你可以得到 ulp2（B）从函数nextafter（B，1.0 / 0.0） - B

Absent overflow and underflow, you can get ulp1(b) from b - nextafter(b, -1.0/0.0). You can get ulp2(b) from nextafter(b, 1.0/0.0) - b.

如果 B / 2'= A＆LT; = 2B ，然后Sterbenz定理告诉你， B - A 是准确的。因此，（B - A） - ULP1 / 2 将是最接近双击来的下限和（b - A）+ ulp2 / 2 将是最接近双击的上限。之前和之后尝试这些价值观和价值观，并挑选工作最宽的区间。

If b/2 <= a <= 2b, then Sterbenz's theorem tells you that b - a is exact. So (b - a) - ulp1 / 2 will be the closest double to the lower bound and (b - a) + ulp2 / 2 will be the closest double to the upper bound. Try these values, and the values immediately before and after, and pick the widest interval that works.

如果 B＆GT; 2A ， B - A＆GT; B / 2 。 B的计算值 - 一个至多一半的ULP关闭。其中 ULP1 是最多两个ULP，作为一个 ulp2 ，所以你给的合理区间是最多两个ULP宽。弄清楚这五个最接近值 B-A 的工作。

If b > 2a, b - a > b/2. The computed value of b - a is off by at most half an ulp. One ulp1 is at most two ulp, as is one ulp2, so the rational interval you gave is at most two ulp wide. Figure out which of the five closest values to b-a work.

如果 A＆GT; 2B ，的ULP B-A 至少是一样大的的ULP b ;如果有什么工作，我敢打赌，它会必须是三个最接近值 B-A 之间。我想其中 A 和 B 有不同的符号同样工作的情况。

If a > 2b, an ulp of b-a is at least as big as an ulp of b; if anything works, I bet it'll have to be be among the three closest values to b-a. I imagine the case where a and b have different signs works similarly.

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

I wrote a small pile of C++ code implementing this idea. It didn't fail random fuzz testing (in a few different ranges) before I got bored of waiting. Here it is:

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;
    }
  }
}

这篇关于最快的算法来识别的最小和最大X中使双precision方程x + A == B真的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！