Find root of a function within an interval adapted from Scipy's bisect.

    Basic bisection routine to find a zero of the function `f` between the
    arguments `a` and `b`. `f(a)` and `f(b)` cannot have the same signs.

    `f` must be jitted via numba.

    Parameters
    ----------
    f : jitted and callable
        Python function returning a number.  `f` must be continuous.
    a : number
        One end of the bracketing interval [a,b].
    b : number
        The other end of the bracketing interval [a,b].
    args : tuple, optional(default=())
        Extra arguments to be used in the function call.
    xtol : number, optional(default=2e-12)
        The computed root ``x0`` will satisfy ``np.allclose(x, x0,
        atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
        parameter must be nonnegative.
    rtol : number, optional(default=4*np.finfo(float).eps)
        The computed root ``x0`` will satisfy ``np.allclose(x, x0,
        atol=xtol, rtol=rtol)``, where ``x`` is the exact root.
    maxiter : number, optional(default=100)
        Maximum number of iterations.
    disp : bool, optional(default=True)
        If True, raise a RuntimeError if the algorithm didn't converge.

    Returns
    -------
    results : namedtuple

def bisect(f, a, b, args=(), xtol=_xtol,
           rtol=_rtol, maxiter=_iter, disp=True):
    """
    Find root of a function within an interval adapted from Scipy's bisect.

    Basic bisection routine to find a zero of the function `f` between the
    arguments `a` and `b`. `f(a)` and `f(b)` cannot have the same signs.

    `f` must be jitted via numba.

    Parameters
    ----------
    f : jitted and callable
        Python function returning a number.  `f` must be continuous.
    a : number
        One end of the bracketing interval [a,b].
    b : number
        The other end of the bracketing interval [a,b].
    args : tuple, optional(default=())
        Extra arguments to be used in the function call.
    xtol : number, optional(default=2e-12)
        The computed root ``x0`` will satisfy ``np.allclose(x, x0,
        atol=xtol, rtol=rtol)``, where ``x`` is the exact root. The
        parameter must be nonnegative.
    rtol : number, optional(default=4*np.finfo(float).eps)
        The computed root ``x0`` will satisfy ``np.allclose(x, x0,
        atol=xtol, rtol=rtol)``, where ``x`` is the exact root.
    maxiter : number, optional(default=100)
        Maximum number of iterations.
    disp : bool, optional(default=True)
        If True, raise a RuntimeError if the algorithm didn't converge.

    Returns
    -------
    results : namedtuple

    """

    if xtol <= 0:
        raise ValueError("xtol is too small (<= 0)")

    if maxiter < 1:
        raise ValueError("maxiter must be greater than 0")

    # Convert to float
    xa = a * 1.0
    xb = b * 1.0

    fa = f(xa, *args)
    fb = f(xb, *args)
    funcalls = 2
    root, status = _bisect_interval(xa, xb, fa, fb)

    # Check for sign error and early termination
    if status == _ECONVERGED:
        itr = 0
    else:
        # Perform bisection
        dm = xb - xa
        for itr in range(maxiter):
            dm *= 0.5
            xm = xa + dm
            fm = f(xm, *args)
            funcalls += 1

            if fm * fa >= 0:
                xa = xm

            if fm == 0 or abs(dm) < xtol + rtol * abs(xm):
                root = xm
                status = _ECONVERGED
                itr += 1
                break

    if disp and status == _ECONVERR:
        raise RuntimeError("Failed to converge")

    return _results((root, funcalls, itr, status))