Calculate the XL-mHG p-value using "Algorithm 2".

    Parameters
    ----------
    N: int
        The length of the list.
    K: int
        The number of 1's in the list.
    X: int
        The XL-mHG ``X`` parameter.
    L: int
        The XL-mHG ``L`` parameter.
    stat: float
        The XL-mHG test statistic.
    tol: float, optional
        The tolerance used for comparing floats. [1e-12]

    Returns
    -------
    float
        The XL-mHG p-value. NaN if floating point precision was insufficient
        for calculating the p-value.

def get_xlmhg_pval2(N, K, X, L, stat, tol=DEFAULT_TOL):
    """Calculate the XL-mHG p-value using "Algorithm 2".

    Parameters
    ----------
    N: int
        The length of the list.
    K: int
        The number of 1's in the list.
    X: int
        The XL-mHG ``X`` parameter.
    L: int
        The XL-mHG ``L`` parameter.
    stat: float
        The XL-mHG test statistic.
    tol: float, optional
        The tolerance used for comparing floats. [1e-12]

    Returns
    -------
    float
        The XL-mHG p-value. NaN if floating point precision was insufficient
        for calculating the p-value.
    """
    # type checking
    assert isinstance(N, int)
    assert isinstance(X, int)
    assert isinstance(L, int)
    assert isinstance(stat, float)
    assert isinstance(tol, float)

    # raise exceptions for invalid parameters
    if not (N >= 1):
        raise ValueError('Invalid value N=%d; must be >= 1.' %(N))
    if not (1 <= X <= N):
        raise ValueError(
            'Invalid value X=%d; must be >= 1 and <= %d.' %(X, N)
        )
    if not (1 <= L <= N):
        raise ValueError(
            'Invalid value L=%d; must be >= 1 and <= %d.' %(L, N)
        )
    if not (0 < stat <= 1.0):
        raise ValueError(
            'Invalid value s=%.1e; must be in (0,1].' %(stat)
        )
    if not (0.0 <= tol < 1.0):
        raise ValueError('Invalid value tol=%.1e; must be in [0,1)' %(tol))

    # special case: stat = 1.0 => pval = 1.0
    if stat == 1.0:
        return 1.0

    W = N-K
    table = np.empty((K+1, L+1), dtype=np.float64)
    table[0,0] = 1.0 # initially, *all* paths have never entered R before
    
    pval = 0.0
    p_start = 1.0
    p = None
    hgp = None
    k = None
    w = None

    # fill dynamic programming table and calculate XL-mHG p-value
    # note: we only need to go over the first L cutoffs, since lower cutoffs
    #       cannot be in R (by definition)
    for n in range(1, L+1):
        
        if K >= n:
            k = n
            p_start *= (float(K-n+1) / float(N-n+1))
        else:
            k = K
            p_start *= (float(n) / float(n-K))

        if p_start == 0.0:
            # not enough floating point precision to calculate
            # the hypergeometric p-value
            return float('nan')

        p = p_start
        hgp = p
        w = n - k
        
        if k == K and (hgp > stat and not is_equal(hgp, stat, tol)):
            # We've exited R (or we were never in it).
            # That means we're done here!
            break
            
        # Check if we've reached R. If so, "go down the diagonal" until we exit R.
        # Three conditions:
        # 1. k >= X         // No configuration with k < X can be in R.
        # 2. w < W          // No configuration with w = W can be in R.
        # 3. pval <= s      // The basic criterion for being in R.
        while k >= X and w < W and (hgp < stat or is_equal(hgp, stat, tol)):
            # We're in R!
            # Note:
            #   For w = W, we always have hgp = 1.0. Since stat < 1.0,
            #   we could just assume that w < W. But this assumption might fail
            #   due to limited floating point accuracy.
            
            # First things first: set table[k, w] to 0 to indicate that this is
            # R territory.
            table[k, w] = 0
            
            # check if we've "just entered" R (this is only possible "from below")
            if table[k-1, w] > 0:
                # calculate the fraction of "fresh" paths (paths which have never entered R before)
                # that enter here, and add that number to r
                pval += (table[k-1, w] * (float(K-k+1)/float(N-n+1)))
                
            p *= (float(k*(N-K-n+k)) / float((n-k+1)*(K-k+1)))
            hgp += p
            w += 1
            k -= 1
            
        # now we're no longer in R
        while k >= 0 and w <= W:
            if k == 0:
                # paths only come in "from the left"
                table[k, w] = table[k, w-1] * (float(W-w+1)/float(N-n+1))
            elif w == 0:
                # paths only come in "from below"
                table[k, w] = table[k-1, w] * (float(K-k+1)/float(N-n+1))
            else:
                # paths come in "from the left" and "from below"
                table[k, w] = table[k, w-1] * (float(W-w+1)/float(N-n+1)) + \
                        table[k-1, w] * (float(K-k+1)/float(N-n+1))
            w += 1
            k -= 1
            
    return pval