This is equivalent to applying the "pset" algorithm once; it's simply expressed in the notation of matrices. And, further:
[1 1]^n x [F_1] = [F_n+1]
[1 0] [F_0] [F_n ]
This is particularly useful, because one can calculate the nth power of a matrix (or of anything) in O(log n) steps with exponentiation by squaring.[0] If we represent a matrix as a list of lists, then the following definitions will work[1][2]...
(def mat-trans (x) ;transpose; this is incredibly terse
(apply map list x))
(def mat-mul (a b)
(let u mat-trans.b
(map (fn (xa)
(map (fn (xb)
(reduce + (map * xa xb)))
u))
a)))
(def mat-id (n) ;n x n identity matrix
(mapn (fn (x)
(mapn (fn (y)
(if (is x y) 1 0))
1 n))
1 n))
(def mat-expt (x n)
(xloop (x x n n total (mat-id (len x)))
(if (is n 0)
total
(even n)
(next (mat-mul x x) (/ n 2) total)
(next x (- n 1) (mat-mul x total)))))
(def fib (n)
(let u (mat-mul (mat-expt '((1 1) (1 0)) n)
'((1) (0)))
u.1.0))
This strategy has the advantage that, once you get it, you can use it for any linear recurrence. For example, if we had
. Also, if one wants, say, these sequences mod m, then it's a fairly simple matter to come up with a "mod"ified version of mat-mul and use that instead--and then the calculation is actually O(log n), because the numbers stay small.
[0] Of course this can be done without the language of matrices. Saying that [1 1; 1 0]^3 = [3 2; 2 1] is the same as saying that the transformation "b -> a+b, a -> b" applied three times is the transformation "b -> 3b+2a, a -> 2b+a". Matrices, and matrix multiplication, are simply a convenient way to represent and compose these transformations.
In fact, as far as I'm concerned, that is the single thing they're good for; and in high school, when I was taught about matrices but not about this stuff, I found it a strange and unmotivated exercise, and so I retained basically none of it, even after it was covered in three consecutive math classes. And I am not known for forgetting or failing to understand things; I've been one of the top students in all of my math classes (except those at an extremely hardcore math camp). I can only imagine how it was for the other students. </rant>
[1] Library functions (with these, the above should work in standard Arc):
(def mapn (f a b)
(if (> a b)
nil
(cons (f a) (mapn f (+ a 1) b))))
(mac xloop (binds . body)
`((rfn next ,(map car pair.binds)
,@body)
,@(map cadr pair.binds)))
[2] It's possible to be slightly more clever than this. E.g. the powers of the [1 1; 1 0] will always look like this:
[F_n+2 F_n+1]
[F_n+1 F_n ]
In theory, you could represent the whole matrix using just two variables (say, F_n+1 and F_n), and the entire computation could be done with four variables.
