What is the computational complexity of Romberg Integration?

Question

I wrote some code, and I'd like to identify the computational complexity. There is one outer loop and two inner loops.

The outer loop goes from i = 1 to some $n$, where $n$ is the number of times a trapezoid is evenly divided. This division results in $2^n + 1$ evaluation points ($k = 1 \dots 2^n + 1$) in the first inner loop.

The second inner loop depends goes from goes from $j = 1 \dots i$, where $i$ is the current index of the outer loop.

So normally, I believe you look at the number of nested loops and that tells you the order of the complexity. However, in this case there are nested loops with very different stopping values. What technique would I use to figure this out, when the stop values are not uniform?

def integrate(func, a, b, steps):
    """
    Uses Romberg Integration to estimate up to best guess R(n, m)

    :param func: The function to integrate
    :param a: start
    :param b: stop
    :param steps: the number of times to divide the interval evenly
    :return: a triangular table of values, where the diagonal contains
    increasing accurate estimates.
    """

    r = np.zeros((steps, steps))
    h = (b - a)
    # Compute Trapezoid Rule estimate:
    r[0, 0] = 0.5*h*(func(a) + func(b))

    # Compute Composite Trap estimates:
    # 2^n is the number of segments
    # N = 2^n + 1 is the number of evaluation points
    # h = (b - a)/2^n is the step size
    for n in range(1, steps):
        h = (b - a) / 2**n

        r[n, 0] = 0.5 * r[n - 1, 0] \
            + h * np.sum(func(a + (2*k - 1)*h)
                         for k in np.arange(1, 2**(n - 1) + 1))

        # Use Richardson extrapolation for higher order accuracy:
        for m in range(1, n + 1):
            r[n, m] = r[n, m - 1] + 1 / (4**m - 1) * (r[n, m - 1] - r[n - 1, m - 1])

    print('The best estimate with convergence order {} and {} parts is {}'.format(steps, 2 ** steps, r[-1, -1]))
    return r