def expm_expl(a: ndarray, output: ndarray) -> tuple[int, int]:
    """
    Implements Taylor expansion, algorithm T_18
    in "Computing the Matrix Exponential with an
    Optimized Taylor Polynomial Approximation",
    Philipp Bader et. al.,
    which minimizes the number of matrix products
    for given number of terms in the expansion.
    """

    tol_m = ExpmConstants.theta_m  # may vary w/ degree, m, in future impls.

    # L1 norm of matrix input:
    l1_norm_a = norm(a, 1)

    requires_scaling = l1_norm_a > tol_m

    s = 0
    A = a
    m = 18

    if requires_scaling:
        s = np.maximum(1, int(np.ceil(np.log2(l1_norm_a / tol_m))))
        #
        # scale `a` by sfactor = 1.0/2^s = 2^(-s):
        #
        sfactor = np.power(2.0, s)
        #
        # A' <- A / sfactor
        #
        A = a / sfactor

    EYE = eye(A.shape[0], dtype=A.dtype)
    A2 = matmul(A, A)
    A3 = matmul(A2, A)
    A6 = matmul(A3, A3)
    B1 = (
        ExpmConstants.a11 * A + ExpmConstants.a21 * A2 + ExpmConstants.a31 * A3
    )
    B2 = (
        ExpmConstants.b11 * A
        + ExpmConstants.b21 * A2
        + ExpmConstants.b31 * A3
        + ExpmConstants.b61 * A6
    )
    B3 = (
        ExpmConstants.b02 * EYE
        + ExpmConstants.b12 * A
        + ExpmConstants.b22 * A2
        + ExpmConstants.b32 * A3
        + ExpmConstants.b62 * A6
    )
    B4 = (
        ExpmConstants.b03 * EYE
        + ExpmConstants.b13 * A
        + ExpmConstants.b23 * A2
        + ExpmConstants.b33 * A3
        + ExpmConstants.b63 * A6
    )
    B5 = (
        ExpmConstants.b24 * A2
        + ExpmConstants.b34 * A3
        + ExpmConstants.b64 * A6
    )

    A9 = B4 + matmul(B1, B5)
    B39 = B3 + A9

    output[:] = B2 + matmul(B39, A9)

    # if scaling by 1/2^s was done then
    # squaring s times is necessary:
    #
    if requires_scaling:
        for j in range(s):
            output[:] = matmul(output, output)

    return (m, s)

    Implements Taylor expansion, algorithm T_18
    in "Computing the Matrix Exponential with an
    Optimized Taylor Polynomial Approximation",
    Philipp Bader et. al.,
    which minimizes the number of matrix products
    for given number of terms in the expansion.
    
cupynumeric/linalg/linalg.py