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.
 
*Tomek Czajka(of topcoder fame) pointed out that my final version was buggy as well. I chose the number of iterations to be 120 but that’s way too small. It doesn't work for c = 1e60.*
 
A double is represented by the mantissa which consists of 52 bits and the exponent which contains 11 bits. One loop iteration either decreases the exponent of our interval by 1 or we find out a new bit of the mantissa. The maximum value of the exponent is 2^11 and the mantissa has 52 bits. Thus we need 2100 steps to figure out the answer.

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:

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

    for iter in range(0, 2100):

    for iter in range(0, 120):

        mid = lo + (hi - lo) / 2
        if (mid * mid * mid < c):
            lo = mid

        return _cubicRoot(c)
==

No more epsilon! But now, because of cases with large exponents, the code runs pretty slow on all cases. An idea is to stop as soon as we don't decrease the lo..hi interval, instead of doing a constant number of iterations. So here’s this faster version:
 
==code(python) |
def _cubicRoot(c):
    lo, hi = 0.0, max(1, c)
    prev_range_size = hi - lo
    while True:
        mid = lo + (hi - lo) / 2
        if (mid * mid * mid < c):
            lo = mid
        else:
            hi = mid
        if hi - lo == prev_range_size:
            break
        prev_range_size = hi - lo
 
    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.

This is still slow when we’re dealing with large numbers. To address that one can binary search on the exponent of the result and then on the mantissa of the result.  Try it out in the comment section.

*BTW Tomek Czajka(of topcoder fame) pointed out that my final version is buggy as well. It doesn't work for c = 1e60. The number of iterations is too small, it needs about 1100 steps.*

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

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

notes:
* Thanks Tomek for the feedback.
* We’ve addressed negative numbers, numbers in (0, 1), overflow problems, absolute and relative error problems.
* Instead of |(a - b) / b| < eps we can use |(a - b) / max(a, b)| < eps. Since a, b >=0 in our problem, we won’t get the nasty cases I was talking about previously, except for the c = 0 case.
* Some tests for your own code -1, 0, 1, 2, 8, 0.001, 1e30, 1e60
* In practice faster methods like 'Newton Rapson method':http://en.wikipedia.org/wiki/Newton's_method  are used.

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