def _inverse_squaring_helper(T0, theta):
    """
    A helper function for inverse scaling and squaring for Pade approximation.

    Parameters
    ----------
    T0 : (N, N) array_like upper triangular
        Matrix involved in inverse scaling and squaring.
    theta : indexable
        The values theta[1] .. theta[7] must be available.
        They represent bounds related to Pade approximation, and they depend
        on the matrix function which is being computed.
        For example, different values of theta are required for
        matrix logarithm than for fractional matrix power.

    Returns
    -------
    R : (N, N) array_like upper triangular
        Composition of zero or more matrix square roots of T0, minus I.
    s : non-negative integer
        Number of square roots taken.
    m : positive integer
        The degree of the Pade approximation.

    Notes
    -----
    This subroutine appears as a chunk of lines within
    a couple of published algorithms; for example it appears
    as lines 4--35 in algorithm (3.1) of [1]_, and
    as lines 3--34 in algorithm (4.1) of [2]_.
    The instances of 'goto line 38' in algorithm (3.1) of [1]_
    probably mean 'goto line 36' and have been intepreted accordingly.

    References
    ----------
    .. [1] Nicholas J. Higham and Lijing Lin (2013)
           "An Improved Schur-Pade Algorithm for Fractional Powers
           of a Matrix and their Frechet Derivatives."

    .. [2] Awad H. Al-Mohy and Nicholas J. Higham (2012)
           "Improved Inverse Scaling and Squaring Algorithms
           for the Matrix Logarithm."
           SIAM Journal on Scientific Computing, 34 (4). C152-C169.
           ISSN 1095-7197

    """
    if len(T0.shape) != 2 or T0.shape[0] != T0.shape[1]:
        raise ValueError('expected an upper triangular square matrix')
    n, n = T0.shape
    T = T0

    # Find s0, the smallest s such that the spectral radius
    # of a certain diagonal matrix is at most theta[7].
    # Note that because theta[7] < 1,
    # this search will not terminate if any diagonal entry of T is zero.
    s0 = 0
    tmp_diag = np.diag(T)
    if np.count_nonzero(tmp_diag) != n:
        raise Exception('internal inconsistency')
    while np.max(np.absolute(tmp_diag - 1)) > theta[7]:
        tmp_diag = np.sqrt(tmp_diag)
        s0 += 1

    # Take matrix square roots of T.
    for i in range(s0):
        T = _sqrtm_triu(T)

    # Flow control in this section is a little odd.
    # This is because I am translating algorithm descriptions
    # which have GOTOs in the publication.
    s = s0
    k = 0
    d2 = _onenormest_m1_power(T, 2) ** (1/2)
    d3 = _onenormest_m1_power(T, 3) ** (1/3)
    a2 = max(d2, d3)
    m = None
    for i in (1, 2):
        if a2 <= theta[i]:
            m = i
            break
    while m is None:
        if s > s0:
            d3 = _onenormest_m1_power(T, 3) ** (1/3)
        d4 = _onenormest_m1_power(T, 4) ** (1/4)
        a3 = max(d3, d4)
        if a3 <= theta[7]:
            j1 = min(i for i in (3, 4, 5, 6, 7) if a3 <= theta[i])
            if j1 <= 6:
                m = j1
                break
            elif a3 / 2 <= theta[5] and k < 2:
                k += 1
                T = _sqrtm_triu(T)
                s += 1
                continue
        d5 = _onenormest_m1_power(T, 5) ** (1/5)
        a4 = max(d4, d5)
        eta = min(a3, a4)
        for i in (6, 7):
            if eta <= theta[i]:
                m = i
                break
        if m is not None:
            break
        T = _sqrtm_triu(T)
        s += 1

    # The subtraction of the identity is redundant here,
    # because the diagonal will be replaced for improved numerical accuracy,
    # but this formulation should help clarify the meaning of R.
    R = T - np.identity(n)

    # Replace the diagonal and first superdiagonal of T0^(1/(2^s)) - I
    # using formulas that have less subtractive cancellation.
    # Skip this step if the principal branch
    # does not exist at T0; this happens when a diagonal entry of T0
    # is negative with imaginary part 0.
    has_principal_branch = all(x.real > 0 or x.imag != 0 for x in np.diag(T0))
    if has_principal_branch:
        for j in range(n):
            a = T0[j, j]
            r = _briggs_helper_function(a, s)
            R[j, j] = r
        p = np.exp2(-s)
        for j in range(n-1):
            l1 = T0[j, j]
            l2 = T0[j+1, j+1]
            t12 = T0[j, j+1]
            f12 = _fractional_power_superdiag_entry(l1, l2, t12, p)
            R[j, j+1] = f12

    # Return the T-I matrix, the number of square roots, and the Pade degree.
    if not np.array_equal(R, np.triu(R)):
        raise Exception('internal inconsistency')
    return R, s, m
dataset/ETHPy150Open scipy/scipy/scipy/linalg/_matfuncs_inv_ssq.py/_inverse_squaring_helper/original