def nanquantile(
    a: ndarray,
    q: float | Iterable[float] | ndarray,
    axis: int | tuple[int, ...] | None = None,
    out: ndarray | None = None,
    overwrite_input: bool = False,
    method: str = "linear",
    keepdims: bool = False,
) -> ndarray:
    """
    Compute the q-th quantile of the data along the specified axis,
    while ignoring nan values.

    Parameters
    ----------
    a : array_like
        Input array or object that can be converted to an array,
        containing nan values to be ignored.
    q : array_like of float
        Quantile or sequence of quantiles to compute, which must be between
        0 and 1 inclusive.
    axis : {int, tuple of int, None}, optional
        Axis or axes along which the quantiles are computed. The default is
        to compute the quantile(s) along a flattened version of the array.
    out : ndarray, optional
        Alternative output array in which to place the result. It must have
        the same shape as the expected output.
    overwrite_input : bool, optional
        If True, then allow the input array `a` to be modified by
        intermediate calculations, to save memory. In this case, the
        contents of the input `a` after this function completes is
        undefined.
    method : str, optional
        This parameter specifies the method to use for estimating the
        quantile.  The options sorted by their R type
        as summarized in the H&F paper [1]_ are:
        1. 'inverted_cdf'
        2. 'averaged_inverted_cdf'
        3. 'closest_observation'
        4. 'interpolated_inverted_cdf'
        5. 'hazen'
        6. 'weibull'
        7. 'linear'  (default)
        8. 'median_unbiased'
        9. 'normal_unbiased'
        The first three methods are discontinuous.  NumPy further defines the
        following discontinuous variations of the default 'linear' (7.) option:
        * 'lower'
        * 'higher',
        * 'midpoint'
        * 'nearest'
    keepdims : bool, optional
        If this is set to True, the axes which are reduced are left in
        the result as dimensions with size one. With this option, the
        result will broadcast correctly against the original array `a`.

    Returns
    -------
    quantile : scalar or ndarray
        If `q` is a single quantile and `axis=None`, then the result
        is a scalar. If multiple quantiles are given, first axis of
        the result corresponds to the quantiles. The other axes are
        the axes that remain after the reduction of `a`. If the input
        contains integers or floats smaller than ``float64``, the output
        data-type is ``float64``. Otherwise, the output data-type is the
        same as that of the input. If `out` is specified, that array is
        returned instead.

    Raises
    ------
    TypeError
        If the type of the input is complex.

    See Also
    --------
    numpy.nanquantile

    Availability
    --------
    Multiple GPUs, Multiple CPUs

    References
    ----------
    .. [1] R. J. Hyndman and Y. Fan,
       "Sample quantiles in statistical packages,"
       The American Statistician, 50(4), pp. 361-365, 1996
    """

    real_axis: int | None
    axes_set: Sequence[int] = ()
    original_shape = a.shape

    if axis is not None and isinstance(axis, Iterable):
        nrm_axis = normalize_axis_tuple(axis, a.ndim)
        if len(axis) == 1:
            real_axis = nrm_axis[0]
            a_rr = a
        else:
            (real_axis, a_rr) = _reshuffle_reshape(a, nrm_axis)
            # What happens with multiple axes and overwrite_input = True ?
            # It seems overwrite_input is reset to False;
            # But `overwrite_input` doesn't matter for the NaN version of this
            # function
            # overwrite_input = False
        axes_set = nrm_axis
    else:
        real_axis = axis
        a_rr = a
        if real_axis is not None:
            axes_set = normalize_axis_tuple(real_axis, a.ndim)
            real_axis = axes_set[0]

    # ndarray of non-NaNs:
    #
    non_nan_counts = asarray(
        count_nonzero(
            logical_not(isnan(a_rr)),
            axis=real_axis,
        )
    )

    # covers both array-like and scalar cases:
    #
    q_arr = np.asarray(q)

    # in the future k-sort (partition)
    # might be faster, for now it uses sort
    # arr = partition(arr, k = floor(nq), axis = real_axis)
    # but that would require a k-sort call for each `q`!
    # too expensive for many `q` values...
    # if no axis given then elements are sorted as a 1D array
    #
    # replace NaN's by dtype.max:
    #
    arr = where(isnan(a_rr), np.finfo(a_rr.dtype).max, a_rr)
    arr.sort(axis=real_axis)

    if arr.dtype.kind == "c":
        raise TypeError("input array cannot be of complex type")

    # return type dependency on arr.dtype:
    #
    # it depends on interpolation method;
    # For discontinuous methods returning either end of the interval within
    # which the quantile falls, or the other; arr.dtype is returned;
    # else, logic below:
    #
    # if is_float(arr_dtype) && (arr.dtype >= dtype('float64')) then
    #    arr.dtype
    # else
    #    dtype('float64')
    #
    # see https://github.com/numpy/numpy/issues/22323
    #
    if method in [
        "inverted_cdf",
        "closest_observation",
        "lower",
        "higher",
        "nearest",
    ]:
        to_dtype = arr.dtype
    else:
        to_dtype = np.dtype("float64")

        # in case dtype("float128") becomes supported:
        #
        # to_dtype = (
        #     arr.dtype
        #     if (arr.dtype == np.dtype("float128"))
        #     else np.dtype("float64")
        # )

    res = nanquantile_impl(
        arr,
        q_arr,
        non_nan_counts,
        real_axis,
        axes_set,
        original_shape,
        _ORDER_FUNCS[method],
        keepdims,
        to_dtype,
        out,
    )

    if out is not None:
        # out = res.astype(out.dtype) -- conversion done inside impl
        return out
    else:
        return res

    Compute the q-th quantile of the data along the specified axis,
    while ignoring nan values.

    Parameters
    ----------
    a : array_like
        Input array or object that can be converted to an array,
        containing nan values to be ignored.
    q : array_like of float
        Quantile or sequence of quantiles to compute, which must be between
        0 and 1 inclusive.
    axis : {int, tuple of int, None}, optional
        Axis or axes along which the quantiles are computed. The default is
        to compute the quantile(s) along a flattened version of the array.
    out : ndarray, optional
        Alternative output array in which to place the result. It must have
        the same shape as the expected output.
    overwrite_input : bool, optional
        If True, then allow the input array `a` to be modified by
        intermediate calculations, to save memory. In this case, the
        contents of the input `a` after this function completes is
        undefined.
    method : str, optional
        This parameter specifies the method to use for estimating the
        quantile.  The options sorted by their R type
        as summarized in the H&F paper [1]_ are:
        1. 'inverted_cdf'
        2. 'averaged_inverted_cdf'
        3. 'closest_observation'
        4. 'interpolated_inverted_cdf'
        5. 'hazen'
        6. 'weibull'
        7. 'linear'  (default)
        8. 'median_unbiased'
        9. 'normal_unbiased'
        The first three methods are discontinuous.  NumPy further defines the
        following discontinuous variations of the default 'linear' (7.) option:
        * 'lower'
        * 'higher',
        * 'midpoint'
        * 'nearest'
    keepdims : bool, optional
        If this is set to True, the axes which are reduced are left in
        the result as dimensions with size one. With this option, the
        result will broadcast correctly against the original array `a`.

    Returns
    -------
    quantile : scalar or ndarray
        If `q` is a single quantile and `axis=None`, then the result
        is a scalar. If multiple quantiles are given, first axis of
        the result corresponds to the quantiles. The other axes are
        the axes that remain after the reduction of `a`. If the input
        contains integers or floats smaller than ``float64``, the output
        data-type is ``float64``. Otherwise, the output data-type is the
        same as that of the input. If `out` is specified, that array is
        returned instead.

    Raises
    ------
    TypeError
        If the type of the input is complex.

    See Also
    --------
    numpy.nanquantile

    Availability
    --------
    Multiple GPUs, Multiple CPUs

    References
    ----------
    .. [1] R. J. Hyndman and Y. Fan,
       "Sample quantiles in statistical packages,"
       The American Statistician, 50(4), pp. 361-365, 1996
    
cupynumeric/_module/stats_order.py