caveat: will fail for large N due to limited precision of floats... works fine on paper though.
[update] caveat 2: assumes a model of computation where ϕ^N can be computed in constant time, which only works on paper.
ϕ_mult = [[1 5] [1 1]]
This is because ϕ(x+y√5)=(x+5y)+(x+y)√5, and
[[1 5] [1 1]] [x y] = [x+5y x+y]
Now, to compute ϕ^N, we just have to compute [x y] = ϕ_mult^N[1 0] (since [1 0] is the representation of 1=1+0√5). Then ϕ^N=x+y√5. This is all exact, no floating point arithmetic!
It's actually just an example of an algebraic number field [1] -- the idea is that "algebraic numbers" (i.e., roots of polynomials) can be equivalently thought of as vectors and matrices. The simplest example, which everyone is familiar with, is the complex numbers, where i=√(-1). In this case, [x y] represents x+y*i, and multiplication by i corresponds to the matrix i_mult=[[0 1] [-1 0]]. Then i_mult^2= -[[1 0] [0 1]], in other words, multiplication by -1!