The blog post author abandons the closed form solution of (φ^n - ψ^n) / (φ - ψ) citing loss of accuracy to demonstrate the pretty neat matrix solution.
The closed form solution is not quite a dead end. The floating point errors can be avoided if one works in the field of number of the form a + b √5 where a and b are rational.
All one needs to is to maintain the numbers a and b separately. The same doubling trick applies for exponentiation.
This is equivalent to the matrix solution. The matrix has eigenvalues φ and ψ, so working with linear combinations of powers of the matrix is the same as working over Q(φ, ψ) = Q(√5). This is not just an equivalence abstractly mathematically, but computationally the solutions are identical and run within a constant factor of each other.
Indeed (φ^n - ψ^n) / (φ - ψ) is what you'd get if you diagonalize the matrix.
However the field you'd be working over is Q[√5], although really we're interested in the subring Z[φ] here. And since φ is (trivially) an eigenvalue of the matrix, let's call it M, this is equivalent to calculating in Z[M] by the Cayley-Hamilton theorem.
Of course some of these constructions might be easier to work with, and certainly something like Z[x]/(x^2-x-1) allows for fairly straightforward calculation, but you're still doing the same multiplications essentially.
And if you're working in a convenient ring you don't really need the ψ^n term, since you can just use that φ^n = F_{n-1} + F_{n} φ.
The closed form solution is not quite a dead end. The floating point errors can be avoided if one works in the field of number of the form a + b √5 where a and b are rational.
All one needs to is to maintain the numbers a and b separately. The same doubling trick applies for exponentiation.