raganwald / homoiconic

Click here to lend your support to: homoiconic and make a donation at www.pledgie.com !

An experiment in publishing code and words about code on a small scale.

This URL has Read+Write access

homoiconic / 2008-12-12 / fibonacci.md
raw normal view history
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 1 A program to compute the nth Fibonacci number
2 ===
3
a4bfc714 » raganwald 2009-09-23 updated resume link 4 Since [I'm looking for a job](http://reginald.braythwayt.com/RegBraithwaiteGH0909_en_US.pdf "Reginald Braithwaite's Resume") and people often like to ask for a [fizzbuzz](http://weblog.raganwald.com/2007/01/dont-overthink-fizzbuzz.html "Don't Overthink FizzBuzz") program to weed out the folks who can't string together a few lines of code, I thought I'd write up a program to compute the nth Fibonacci number. There's an intriguing bit of matrix math involved, so I learned something while implementing it.
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 5
899d451c » Reginald Braithwaite 2008-12-13 added an unspace credit 6 Recently I was having lunch with [some of Toronto's most interesting Ruby developers](http://unspace.ca), and the subject of interview questions came up. Specifically, writing a program to compute the nth Fibonacci number. Unsurprisingly, we agreed it might be useful as a screener for weeding out the people who were completely delusional about their prospects as a professional programmer. You know, the type of person who stares blankly at the screen and has no idea where to start, even when told that all we wanted was a program that correctly prints an answer. No tests, no specs, no shouldas, no passengers, no CSS, just prove you actually can write something, anything.
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 7
1a3e1d75 » raganwald 2008-12-12 clarified a sentance 8 We also agreed that any [Guess the answer I'm thinking of](http://www.nomachetejuggling.com/2008/12/11/my-least-favorite-interview-question/ "My Least Favorite Interview Question » Absolutely No Machete Juggling") problem is terrible, including Fibonacci. Should it iterate? Recurse? Memoize? Is it overkill to use advanced math to be really fast? Or will you get tossed out of the interview for writing a naive recursive solution? Unless the interviewer is very specific that they just want you to prove that you actually have written ten or more lines of working Ruby code in your life, there could be any number of reasons to disqualify any solution.
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 9
10 But we batted around a few ideas for how to gild the lily, and since I have just found myself contemplating the prospect of being asked to write a program to compute the nth Fibonacci number, I thought I'd get some practice and write one in advance. Since I'm doing it for myself, I'm going to optimize for maximum personal growth.
11
12 Here we go...
13
14 Enter The Matrix
15 ---
16
e6c44800 » Reginald Braithwaite 2008-12-13 fixed some inaccuracies abo... 17 One problem with calculating a Finonacci number is that naive algorithms require _n_ addition operations. There are some interesting things we can do to improve on this. One way is to transform _n_ additions into raising something to the power of *n*. This turns _n_ additions into _n_ multiplications. That seems retrograde, but hold on to your disbelief.
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 18
e6c44800 » Reginald Braithwaite 2008-12-13 fixed some inaccuracies abo... 19 This is actually nice, because there is a trick about raising something to a power that we can exploit. But first things first. As explained in [Sum of Fibonacci numbers?](http://expertvoices.nsdl.org/cornell-cs322/2008/03/25/sum-of-fibonacci-numbers/), we can express the Fibonacci number `F(n)` using a 2x2 matrix:
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 20
21 [ 1 1 ] n [ F(n+1) F(n) ]
22 [ 1 0 ] = [ F(n) F(n-1) ]
23
a94dd2da » raganwald 2008-12-14 Matrix Math correct ion fro... 24 [Multiplying two matrices](http://www.maths.surrey.ac.uk/explore/emmaspages/option1.html "Matrices and Determinants") is a little interesting if you have never seen it before:
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 25
a94dd2da » raganwald 2008-12-14 Matrix Math correct ion fro... 26 [ a b ] [ e f ] [ ae + bg af + bh ]
27 [ c d ] times [ g h ] = [ ce + dg cf + dh ]
28
29 Our matrices always have diagonal symmetry, so we can simplify the calculation because _c_ is always equal to _b_:
30
31 [ a b ] [ d e ] [ ad + be ae + bf ]
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 32 [ b c ] times [ e f ] = [ bd + ce be + cf ]
33
ee1fac1f » raganwald 2008-12-14 proved the commutative prop... 34 Now we are given that we are multiplying two matrices with diagonal symmetry. Will the result have diagonal symmetry? In other words, will `ae + bf` always be equal to `bd + ce`? Remember that `a = b + c` at all times and `d = e + f` provided that each is a power of `[[1,1], [1,0]]`. Therefore:
35
36 ae + bf = (b + c)e + bf
37 = be + ce + bf
38 bd + ce = b(e + f) + ce
39 = be + bf + ce
40
41 That simplifies things for us, we can say:
42
33808d75 » raganwald 2008-12-14 fixed a place where some co... 43 [ a b ] [ d e ] [ ad + be ae + bf ]
44 [ b c ] times [ e f ] = [ ae + bf be + cf ]
ee1fac1f » raganwald 2008-12-14 proved the commutative prop... 45
46 And thus, we can always work with three elements instead of four. Let's express this as operations on arrays:
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 47
a94dd2da » raganwald 2008-12-14 Matrix Math correct ion fro... 48 [a, b, c] times [d, e, f] = [ad + be, ae + bf, be + cf]
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 49
ee1fac1f » raganwald 2008-12-14 proved the commutative prop... 50 Which we can code in Ruby:
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 51
52 def times(*ems)
53 ems.inject do |product, matrix|
54 a,b,c = product; d,e,f = matrix
a94dd2da » raganwald 2008-12-14 Matrix Math correct ion fro... 55 [a*d + b*e, a*e + b*f, b*e + c*f]
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 56 end
57 end
58
59 times([1,1,0]) # => [1, 1, 0]
60 times([1,1,0], [1,1,0]) # => [2, 1, 1]
61 times([1,1,0], [1,1,0], [1,1,0]) # => [3, 2, 1]
62 times([1,1,0], [1,1,0], [1,1,0], [1,1,0]) # => [5, 3, 2]
63 times([1,1,0], [1,1,0], [1,1,0], [1,1,0], [1,1,0]) # => [8, 5, 3]
64
65 Very interesting. We could write out a naive implementation that constructs a long array of copies of `[1,1,0]` and then calls `times`:
66
67 def naive_power(m, n)
68 times(*(1..n).map { |n| m })
69 end
70
71 naive_power([1,1,0], 1) # => [1, 1, 0]
72 naive_power([1,1,0], 2) # => [2, 1, 1]
73 naive_power([1,1,0], 3) # => [3, 2, 1]
74 naive_power([1,1,0], 4) # => [5, 3, 2]
75 naive_power([1,1,0], 5) # => [8, 5, 3]
76
49ae2542 » Reginald Braithwaite 2008-12-12 a few fixes 77 Now let's make an observation: instead of accumulating a product by iterating over the list, let's [Divide and Conquer](http://www.cs.berkeley.edu/~vazirani/algorithms/chap2.pdf). Let's take the easy case: Don't you agree that `times([1,1,0], [1,1,0], [1,1,0], [1,1,0])` is equal to `times(times([1,1,0], [1,1,0]), times([1,1,0], [1,1,0]))`? And that this saves us an operation, since `times([1,1,0], [1,1,0], [1,1,0], [1,1,0])` is implemented as:
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 78
79 times([1,1,0],
80 times([1,1,0],
81 times([1,1,0],[1,1,0]))
82
83 Whereas `times(times([1,1,0], [1,1,0]), times([1,1,0], [1,1,0]))` can be implemented as:
84
85 begin
86 double = times([1,1,0], [1,1,0])
87 times(double, double)
88 end
89
90 This only requires two operations rather than three. Furthermore, it is recursive. `naive_power([1,1,0], 8)` requires seven operations. However, it can be formulated as:
91
92 begin
93 quadruple = begin
94 double = times([1,1,0], [1,1,0])
95 times(double, double)
96 end
97 times(quadruple, quadruple)
98 end
99
100 Now we only need three operations compared to seven. Of course, we left out how to deal with odd numbers. Fixing that also fixes how to deal with even numbers that aren't neat powers of two:
101
102 def power(m, n)
103 if n == 1
104 m
105 else
106 halves = power(m, n / 2)
107 if n % 2 == 0
108 times(halves, halves)
109 else
110 times(halves, halves, m)
111 end
112 end
113 end
114
115 power([1,1,0], 1) # => [1, 1, 0]
116 power([1,1,0], 2) # => [2, 1, 1]
117 power([1,1,0], 3) # => [3, 2, 1]
118 power([1,1,0], 4) # => [5, 3, 2]
119 power([1,1,0], 5) # => [8, 5, 3]
120
121 Now we can write our complete fibonacci function:
122
123 def fib(n)
d85aab9c » raganwald 2008-12-12 fixed a bug 124 return n if n < 2
125 power([1,1,0], n - 1).first
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 126 end
127
128 And dress things up in idiomatic Ruby using the anonymous module pattern:
129
130 class Integer
131
132 include(Module.new do
133
134 times = lambda do |*ems|
135 ems.inject do |product, matrix|
136 a,b,c = product; d,e,f = matrix
a94dd2da » raganwald 2008-12-14 Matrix Math correct ion fro... 137 [a*d + b*e, a*e + b*f, b*e + c*f]
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 138 end
139 end
140
141 power = lambda do |m, n|
142 if n == 1
143 m
144 else
145 halves = power.call(m, n / 2)
146 if n % 2 == 0
147 times.call(halves, halves)
148 else
149 times.call(halves, halves, m)
150 end
151 end
152 end
153
154 define_method :matrix_fib do
155 return self if self < 2
156 power.call([1,1,0], self - 1).first
157 end
158
159 end)
160
161 end
162
49ae2542 » Reginald Braithwaite 2008-12-12 a few fixes 163 (0..20).map { |n| n.matrix_fib }
164 # => [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765]
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 165
166 We're done!
167
94aa33b7 » raganwald 2009-02-02 Updated links, especially t... 168 p.s. No we're not done: [Another program to compute the nth Fibonacci number](http://github.com/raganwald/homoiconic/tree/master/2008-12-17/another_fibonacci.md#readme).
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 169
6c25f36b » raganwald 2008-12-17 linked the first fibonacci ... 170 p.p.s. No, this isn't [the fastest implementation](http://blade.nagaokaut.ac.jp/cgi-bin/scat.rb/ruby/ruby-talk/194815 "Fast Fibonacci method") by far. But it beats the pants off of a naive iterative implementation. See [fibonacci.rb](http:fibonacci.rb) for details.
c017b61d » Reginald Braithwaite 2008-12-13 pulled comments from Christ... 171
6c25f36b » raganwald 2008-12-17 linked the first fibonacci ... 172 ---
173
174 Christoph Temmel [forked homoiconic](http://github.com/kar8nga/homoiconic/tree "kar8nga's homoiconic at master &mdash; GitHub") and added these observations:
175
176 *Another option would be to decompose *
88cb589b » Christoph Temmel 2008-12-13 comment on matrices 177
178 A = [ 1 1 ]
179 [ 1 0 ]
180
6c25f36b » raganwald 2008-12-17 linked the first fibonacci ... 181 *into along its eigenspace (eigenvalues \lambda\_{1,2}=-\frac{1\pm\sqrt{5}}{2} = -\frac{1}{2}\mp\frac{\sqrt{5}{2}}) to get*
88cb589b » Christoph Temmel 2008-12-13 comment on matrices 182
72b150a8 » Christoph Temmel 2008-12-13 comment on matrices 183 A = Q^t D
88cb589b » Christoph Temmel 2008-12-13 comment on matrices 184
6c25f36b » raganwald 2008-12-17 linked the first fibonacci ... 185 *and*
72b150a8 » Christoph Temmel 2008-12-13 comment on matrices 186
187 D = [ \lambda_1 0 ]
188 [ 0 \lambda_2 ]
189
6c25f36b » raganwald 2008-12-17 linked the first fibonacci ... 190 *where Q is the orthonormal matrix of the normated eigenvectors with QQ^t = I and D is the above diagonal matrix with the eigenvalues \lambda\_{1,2} (see also [Wikipedia](http://en.wikipedia.org/wiki/Symmetric_matrix#Properties)). Now it's easy to take powers of A*
88cb589b » Christoph Temmel 2008-12-13 comment on matrices 191
192 A^n = Q^t D^n D
193
6c25f36b » raganwald 2008-12-17 linked the first fibonacci ... 194 *The only difficulty is, that in order to avoid floating point arithmetic one would have to do the powers in the field \mathbb{Q}[\sqrt{5}] - this would ask for a custom datatype with overloading of addition and multiplication (the part needed for this exercise). *
88cb589b » Christoph Temmel 2008-12-13 comment on matrices 195
196
197 ----
0e6f843e » Reginald Braithwaite 2008-12-12 A program to compute the nt... 198
199 Subscribe to [new posts and daily links](http://feeds.feedburner.com/raganwald "raganwald's rss feed"): <a href="http://feeds.feedburner.com/raganwald"><img src="http://feeds.feedburner.com/~fc/raganwald?bg=&amp;fg=&amp;anim=" height="26" width="88" style="border:0" alt="" align="top"/></a>
200
a4bfc714 » raganwald 2009-09-23 updated resume link 201 Reg Braithwaite: [Home Page](http://reginald.braythwayt.com), [CV](http://reginald.braythwayt.com/RegBraithwaiteGH0909_en_US.pdf ""), [Twitter](http://twitter.com/raganwald)