Numpy: Raise diagonalizable square matrix to infinite power

Question

Consider a Markovian process with a diagonalizable transition matrix A such that A=PDP^-1, where D is a diagonal matrix with eigenvalues of A, and P is a matrix whose columns are eigenvectors of A.

To compute, for each state, the likelihood of ending up in each absorbing state, I'd like to raise the transition matrix to the power of n, with n approaching infinity:

A^n=PD^nP^-1

Which is the pythonic way of doing this in Numpy? I could naively compute the eigenvalues and eigenvectors of A, raising the eigenvalues to infinity. Due to my assumption that I have a transition matrix, we would have only eigenvalues equal to one (which will remain one), and eigenvalues between 0 and 1, which will become zero (inspired by this answer):

import numpy as np
from scipy import linalg    

# compute left eigenvalues and left eigenvectors
eigenvalues, leftEigenvectors = linalg.eig(transitionMatrix, right=False, left=True)

# for stationary distribution, eigenvalues and vectors are real (up to numerical precision)
eigenvalues = eigenvalues.real
leftEigenvectors = leftEigenvectors.real

# create a filter to collect the eigenvalues close to one
absoluteTolerance = 1e-10
mask = abs(eigenvalues - 1) < absoluteTolerance

# raise eigenvalues to the power of infinity
eigenvalues[mask] = 1
eigenvalues[~mask] = 0

D_raised = np.diag(eigenvalues)

A_raised = leftEigenvectors @ D_raised @ linalg.inv(leftEigenvectors)

Would this be the recommended approach, i.e., one that is both numerically stable and efficient?