import dask.array as da

def _triu_indices(n, k=0, m=None):
    # https://stackoverflow.com/questions/242711/algorithm-for-index-numbers-of-triangular-matrix-coefficients

    if m is None:
        m = n

    if k < 0 or m != n:
        # If k is negative or M does not match n then use normal method,
        # could not figure out how to calculate those algebraically:
        return da.nonzero(~da.tri(n, m, k=k - 1, dtype=bool))

    def _triangular_number(n, k=0):
        """Get the triangular number."""
        if k >= 0:
            T_n = (n - k) * (n - k + 1) * 0.5
        else:
            # When k is negative add triangular numbers from the other
            # side instead:
            k = -1 * k + 1
            T_n = (n - 0) * (n - 0 + 1) * 0.5
            T_n += (n - 1) * (n - 1 + 1) * 0.5
            T_n -= (n - k) * (n - k + 1) * 0.5

        return T_n

    def _triangular_root(T_n, k=0):
        """Solve _triangular_number for n instead."""
        n = (da.sqrt(8*T_n + (1 - 2*k)*(1 - 2*k) + 4*k - 4*k*k) - (1 - 2*k)) * 0.5
        return n

    # Number of triangular points:
    T_n = int(_triangular_number(n=n, k=k))
    T_n_range = da.arange(T_n, dtype=int)

    #
    ii = T_n - 1 - T_n_range

    # Triangular root:
    K = da.floor(_triangular_root(ii, k=k))

    #
    jj = ii - _triangular_number(n=K, k=k)

    #
    i = n - 1 - K
    j = m - 1 - jj

    return (i, j)