Our goal is just to find a vector that satisfies ~Pv = v~ and sums to ~1~. By rearranging, we can find that:

$$\displaystyle Pv = v \implies (P-I)v = 0$$

Thus, we just need to find the kernel of ~P-I~. In particular, since the answer exists and is unique, we know that ~\dim \ker(P-I) = 1~ so the basis has exactly one element, and we can just scale the single basis vector so that its sum is ~1~.

This can be done with Gaussian elimination.

Alternatively, we can just solve directly for the linear system ~Pv = 0~ by further constraining it so that the only result is where ~\sum_{i=1}^N v_i = 1~. This can be done by adding an additional row of ~1~s to ~P~, and making the target vector ~\langle 0, 0, \dots, 0, 1 \rangle~.

Time Complexity: ~\mathcal O(N^3)~