@ft.lru_cache
def f(z: complex, c: complex) -> complex:
    """general function f, such that f_c(z) = f(z, c)"""
    # Early exit if z is too large to avoid overflow
    if np.isnan(z) or np.abs(z) > 1e10:
        return np.nan + np.nan * 1j
    try:
        return z**2 + c
    except OverflowError:
        return np.nan + np.nan * 1j

@ft.lru_cache
def v(k: int, c: complex) -> complex:
    r"""
    general recursive function v, such that
    v_c(k) = v(k, c) = f_c(v(k-1, c)) = f(v(k-1, c), c) for all k >= 1, and
    v_c(0) = f_c(0) = f(0, c).
    """
    assert k >= 0, "input k must be non-negative."
    # Early exit if previous value is nan (overflow detected)
    prev = 0 if k == 0 else v(k-1, c)
    if np.isnan(prev) or (isinstance(prev, complex) and np.isnan(prev.real)):
        return np.nan + np.nan * 1j
    return f(0, c) if k == 0 else f(prev, c)


def approximator(c: complex, max_k: int = 50) -> complex:
    """use a max_k to approximate the final bounded value of the sequence {v_c(k)}"""
    return np.absolute(v(max_k, c))


def approx_binner(
    a: np.ndarray,
    pre_apply: callable = np.log2,
    bins: Iterable = range(-2, 308)
):
    """
    Helper function that clusters values of array `a` into different bins and
    returns the indices of the bins to which each value in input array belongs.

    Parameters
    ----------
    pre_apply : callable
        a function applied to input array before binning its values.
    bins : Iterable
        a 1-D list-like object that determins cutoff-point of bins.
    """
    # Mask non-positive values to avoid log2(0) or log2(negative)
    a_masked = np.where(a > 0, a, np.nan)
    return np.digitize(pre_apply(a_masked), bins=bins)

def binary_indicator(c: complex, max_k: int = 50, threshold: float = 100) -> bool:
    r"""
    Indicator function that indicate if the series {v_c(k)} stays bounded, ie,
    belongs to Mandelbrot set (returns True) or not (returns False).
    
    There are rooms for better algorithm implementation since the current
    approach is only a simple approximation to determine is a series is bounded
    but not always correct compared to theory.

    Parameters
    ----------
        c (complex): a fixed complex number for the function v (and thus f).
        max_k (int): the max allowed k to be used to compute the series {v(k, c)}.
        threshold (float, >=0): 
          the customized bound ($\gamma_c$) to determine if the series {v(k, c)}
          where k=0, 1,... max_k, is bounded or not.
    """
    return approximator(c, max_k) <= threshold