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