code: def _fast_numerical_inverse(dist, tol=1e-12, max_intervals=100000):
      """
      Generate fast, approximate PPF (inverse CDF) of probability distribution.
  
      `_fast_numerical_inverse` accepts `dist`, an object representing the
      distribution for which a fast approximate PPF is desired, and returns an
      object `fni` with methods that approximate `dist.ppf` and `dist.rvs`.
      For some distributions, these methods may be faster than those of `dist`
      itself.
  
      Parameters
      ----------
      dist : object
          Object representing distribution for which fast approximate PPF is
          desired; e.g., a frozen instance of `scipy.stats.rv_continuous`.
      tol : float, optional
          u-error tolerance. The default is 1e-12.
      max_intervals : int, optional
          Maximum number of intervals in the cubic Hermite Spline used to
          approximate the percent point function. The default is 100000.
  
      Returns
      -------
      H : scipy.interpolate.CubicHermiteSpline
          Interpolant of the distributions's PPF.
      intervals : int
          The number of intervals of the interpolant.
      midpoint_error : float
          The maximum u-error at an interpolant interval midpoint.
      a, b : float
          The left and right endpoints of the valid domain of the interpolant.
  
      """
      dist, tol, max_intervals = _fni_input_validation(
          dist, tol, max_intervals
      )
  
      # [1] Section 2.1: "For distributions with unbounded domain, we have to
      # chop off its tails at [a] and [b] such that F(a) and 1-F(b) are small
      # compared to the maximal tolerated approximation error."
      p = np.array(
          [dist.ppf(tol / 10), dist.isf(tol / 10)]
      )  # initial interval
  
      # [1] Section 2.3: "We then halve this interval recursively until
      # |u[i+1]-u[i]| is smaller than some threshold value, for example, 0.05."
      u = dist.cdf(p)
      while p.size - 1 <= np.ceil(max_intervals / 2):
          i = np.nonzero(np.diff(u) > 0.05)[0]
          if not i.size:
              break
  
          p_mid = (p[i] + p[i + 1]) / 2
          # Compute only the new values and insert them in the right places
          # [1] uses a linked list; we can't do that efficiently
          u_mid = dist.cdf(p_mid)
          p = np.concatenate((p, p_mid))
          u = np.concatenate((u, u_mid))
          i_sort = np.argsort(p)
          p = p[i_sort]
          u = u[i_sort]
  
      # [1] Section 2.3: "Now we continue with checking the error estimate in
      # each of the intervals and continue with splitting them until [it] is
      # smaller than a given error bound."
      u = dist.cdf(p)
      f = dist.pdf(p)
      while p.size - 1 <= max_intervals:
          # [1] Equation 4-8
          try:
              H = CubicHermiteSpline(u, p, 1 / f)
          except ValueError:
              message = (
                  "The interpolating spline could not be created. This "
                  "is often caused by inaccurate CDF evaluation in a "
                  "tail of the distribution. Increasing `tol` can "
                  "resolve this error at the expense of lower accuracy."
              )
              raise ValueError(message)
          # To improve performance, add update feature to CubicHermiteSpline
  
          # [1] Equation 12
          u_mid = (u[:-1] + u[1:]) / 2
          eu = np.abs(dist.cdf(H(u_mid)) - u_mid)
  
          i = np.nonzero(eu > tol)[0]
          if not i.size:
              break
  
          p_mid = (p[i] + p[i + 1]) / 2
          u_mid = dist.cdf(p_mid)
          f_mid = dist.pdf(p_mid)
          p = np.concatenate((p, p_mid))
          u = np.concatenate((u, u_mid))
          f = np.concatenate((f, f_mid))
          i_sort = np.argsort(p)
          p = p[i_sort]
          u = u[i_sort]
          f = f[i_sort]
  
      # todo: add test for monotonicity [1] Section 2.4
      # todo: deal with vanishing density [1] Section 2.5
      return H, eu, p.size - 1, u[0], u[-1]
  
file_path: catboost/contrib/python/scipy/py3/scipy/stats/_rvs_sampling.py
imports: import numpy
  from scipy._lib._util import check_random_state
  from scipy.interpolate import CubicHermiteSpline
  from scipy import stats
label_desc: Use image node_id 3 for calling a global function with example usage: _fast_numerical_inverse(dist, tol, max_intervals) and returns: H,  eu,  unknown,  u,  u