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

    for iter in range(0, 120):

    for iter in range(0, 1100):

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

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

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?