Calculate the XL-mHG test statistic using recurrence relations.
    
    Parameters
    ----------
    v: np.ndarray with dtype=np.uint8
        The ranked list.
    X: int
        The XL-mHG ``X`` parameter.
    L: int
        The XL-mHG ``L`` parameter.
    tol: float, optional
        The tolerance used for comparing floats. [1e-12]
    
    Returns
    -------
    s: float
        The XL-mHG test statistic.
    n_star: int
        The (first) cutoff at which s was attained.

def get_xlmhg_stat(v, X, L, tol=DEFAULT_TOL):
    """Calculate the XL-mHG test statistic using recurrence relations.
    
    Parameters
    ----------
    v: np.ndarray with dtype=np.uint8
        The ranked list.
    X: int
        The XL-mHG ``X`` parameter.
    L: int
        The XL-mHG ``L`` parameter.
    tol: float, optional
        The tolerance used for comparing floats. [1e-12]
    
    Returns
    -------
    s: float
        The XL-mHG test statistic.
    n_star: int
        The (first) cutoff at which s was attained.
    """
    assert isinstance(v, np.ndarray) and v.ndim == 1 and \
           np.issubdtype(v.dtype, np.integer)
    assert isinstance(X, int)
    assert isinstance(L, int)
    assert isinstance(tol, float)

    N = v.size
    if not N > 0:
        raise ValueError('List is empty!')
    if not (1 <= X <= N):
        raise ValueError(
            'Invalid value X=%d; should be >= 1 and <= %d.' % (X, N)
        )
    if not (1 <= L <= N):
        raise ValueError(
            'Invalid value L=%d; should be >= 1 and <= %d.' % (L, N)
        )
    if not (0.0 <= tol < 1.0):
        raise ValueError('Invalid value tol=%.1e; should be in [0,1)' % (tol))
    
    K = int(np.sum(v != 0))
    if K == 0:
        return 1.0, 0
 
    p = 1.0
    stat = 1.1
    n_star = 0
    k = 0
    for n in range(L):
        if v[n] == 0:
            # calculate f(k; N,K,n+1) from f(k; N,K,n)
            p *= (float((n+1)*(N-K-n+k)) / float((N-n)*(n-k+1)))
        else:
            # we hit a 1
            # calculate f(k+1; N,K,n+1) from f(k; N,K,n)
            p *= (float((n+1)*(K-k)) / float((N-n)*(k+1)))
            k += 1
            # calculate hypergeometric p-value only if enough elements have
            # been seen
            if k >= X:
                hgp = get_hgp(p, k, N, K, n+1)
                if hgp < stat and not is_equal(hgp, stat, tol):
                    stat = hgp
                    n_star = n+1

    stat = min(stat, 1.0)
    return stat, n_star