Fastest algorithm to identify the smallest and largest x that make the double-precision equation x + a == b true

Question

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 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.

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.

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.

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.

Could you use binary search to bisect to the values that you want? It would seem like this should be possible since the number of bits is low. — templatetypedef, Jun 14 '14 at 18:47
@templatetypedef the “128 floating-point additions” solution I sketch out is a binary search over the representation of floating-point numbers, and the one that I don't want to show because I don't know if it is actually an improvement reduces the initial interval to bisect by computing an over-approximated range of candidates, that would then need to be refined by binary search. — Pascal Cuoq, Jun 14 '14 at 18:54
@templatetypedef I am sort of hoping that someone will come up with a theorem of floating-point arithmetics that solves the problem more elegantly. — Pascal Cuoq, Jun 14 '14 at 18:56
Have you considered hybrid approaches? Do some analysis of the numbers to get the approximate range, then refine using binary search. — Patricia Shanahan, Jun 14 '14 at 19:08
@PatriciaShanahan The solution I didn't want to show is to start from the range rd(nextafter(b, -INFINITY) - a) … ru(nextafter(b, +INFINITY) - a) as initial approximation of the lower bound, then use either bitwise operations on the floating-point representation or floating-point arithmetics to refine the result. Then refine the upper bound starting from the same initial interval. — Pascal Cuoq, Jun 14 '14 at 19:18
I'll try to think of something else, but that was my best initial idea. — Patricia Shanahan, Jun 14 '14 at 19:24
@PascalCuoq: Seems like you can get `ulp1` and `ulp2` using a couple `nextafter`s and subtractions. If `a` and `b` have the same sign, Sterbenz's theorem tells you that `b-a` is exact and it's also at most `b`. What about scaling `ulp1` and `ulp2` up by `1 + 0x1p-52` before pushing the bounds out, then just testing two or three nearby values? (Don't have time to flesh this out right now. I guess I'll post it as an answer now and work on it later.) — tmyklebu, Jun 14 '14 at 19:24
@tmyklebu You lost me at “If a and b have the same sign, Sterbenz's theorem tells you that b-a is exact”. What if `a` is 1 and `b` is `0x1.0p80`? The Sterbenz theorem I know is for b/2 ≤ a ≤ 2*b. — Pascal Cuoq, Jun 14 '14 at 19:29

tmyklebu · Accepted Answer · 2014-06-15T01:41:01.993

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

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.

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.

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.

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.

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.

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.

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

Great! Exactly what I was hoping that someone would find. One question though: reading from my phone, I don't see the place where you had to worry about the representation of the empty set when the empty set is indeed the best answer (`x + 1.0 == 0x1.0p-80`, for instance) — Pascal Cuoq, Jun 15 '14 at 09:17
@PascalCuoq: I'm inconsistent about that. Dealing with the NaN/infinity cases, I just return. Later on, I return with `xlo > xhi`. — tmyklebu, Jun 15 '14 at 12:54

Fastest algorithm to identify the smallest and largest x that make the double-precision equation x + a == b true

1 Answers1

Linked