Why your bisection search is wrong

Cosmin Negruseri
19 mai 2014

What is bisection search? The bisection method or bisection search is a numerical algorithm for finding a value x such that f(x) = 0, for a given continuous function f. It works by repeatedly bisecting an interval and choosing a subinterval that contains x. It's pretty simple and robust, but it has few gotchas.

Let's solve the following problem:

For a given number c find it's cubic root using the +, -, *, / operations.

Try solving the problem on your own, before reading below.

Let's choose f(x) = x³ - c. f is continuous and x is the cubic root of c, when f(x) = 0. Thus, we can apply the bisection method.

def cubicRoot(c): {
  lo, hi = 0.0, c
  while lo * lo * lo != c:
    mid = (lo + hi) / 2;
    if mid * mid * mid < c:
       lo = mid
    else:
       hi = mid

  return lo

Any bugs? Well, quite a few. Try to spot as many as you can, before reading on.

You may notice the precision issue right from the start. We'll discuss it a bit later.

What else? The code doesn’t work for negative values of c. This is easily fixable:

def cubicRoot(c): {
  lo, hi = 0.0, c
  while lo * lo * lo != c:
    mid = (lo + hi) / 2;
    if mid * mid * mid < c:
       lo = mid
    else:
       hi = mid

  return lo

def cubicRoot(c):
   if c < 0:
       return - _cubicRoot(-c)
   else:
       return _cubicRoot(c)
}

What else? This code doesn't work for c = 1/1000. Why is that? We’re setting the cubic root upper limit to c (line 3). But, the cubic root of c in (0, 1) is larger than c, meaning the upper bound we’re setting is wrong.

Let's fix it:

def _cubicRoot(c):
    lo, hi = 0.0, max(1, c)

    while lo * lo * lo != c:
        mid = (lo + hi) / 2
        if (mid * mid * mid < c):
            lo = mid
        else:
            hi = mid

    return lo

def cubicRoot(c):
    if c < 0:
        return -_cubicRoot(-c)
    else:
        return _cubicRootABS(c)

Going back back to the precision issue:

If you've read nearly all binary searches and merge sorts are wrong you know that mid = (lo + hi) / 2 has an overflow problem. So we change that line to mid = lo + (hi - lo) / 2.

Testing for equality doesn't work with floating point numbers. The first idea that comes to mind is to use an absolute error comparison (|a - b| < eps).

Instead of:

while lo * lo * lo != c

switch to

while abs(c - lo * lo * lo) < 1e-3:

This doesn't work. For large numbers like 1e37 the code goes into an infinite loop. Try to figure out why. Discuss it in the comment section. Let’s try using the relative error (|(a - b) / b| < eps). There are some weird cases when a and b are close or equal to 0. Can the code be cleaner?

After each while loop iteration we learn something new about x’s range. A double has only 64 bits of precision. So instead of a tricky floating point stopping criteria we can run the loop a fixed number of times so that the interval is as small as the precision of our floating point numbers allows. Now the algorithm is:

def _cubicRoot(c):
    lo, hi = 0.0, max(1, c)

    for iter in range(0, 120):
        mid = lo + (hi - lo) / 2
        if (mid * mid * mid < c):
            lo = mid
        else:
            hi = mid

    return lo

def cubicRoot(c):
    if c < 0:
        return -_cubicRoot(-c)
    else:
        return _cubicRoot(c)

Done! No more epsilon! The result is as accurate as possible given a double representation.

BTW Tomek Czajka pointed out that my final version is buggy as well. It doesn't work for c = 1e60. The problem is the number of iterations is too small, and for it to work correctly it needs about 1100 steps.

We’ve addressed negative numbers, numbers in (0, 1), overflow problems, absolute and relative error problems.

I'm curious, did your solution have any of these problems?