docs: 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))