The title is hyperbole (I'm sure that I've written much, much worse algorithms, many times), but the breakdown of Fibonacci sequence algorithms is really enjoyable.
Typical way to gain views and make it look like you know what you are talking about (although the author does seem to know good algorithm development). Hell, comparing it to Bogosort is a stretch. Bogosort is not even a naive algorithm.
FWIW, in my Programming Languages Theory class the first sorting algorithm we learned for Prolog was Permutation sort, which is a better version of Bogosort. Instead of trying a random permutation each time, Permutation sort will try each permutation once. In Prolog, this is a (the?) naive sort.
The code below consists of declaring that S is a sorted version of L as long as S is a permutation of L and S is sorted.
permutation_sort(L,S) :- permutation(L,S), sorted(S).
sorted([]).
sorted([_]).
sorted([X,Y|ZS]) :- X =< Y, sorted([Y|ZS]).
permutation([],[]).
permutation([X|XS],YS) :- permutation(XS,ZS),select(X,YS,ZS).
Where is ZS defined? Bear in mind that I only took half a module's worth (6 modules in a year at my university) of Prolog and we really just used it for learning predicate logic.
Nevertheless, using Prolog (especially when my code worked!) blew my mind. In the assignment I had to write a program that mimicked part of an aircraft controller in that each aircraft in it's flightplan had to be separated from all other aircraft. Do you think Prolog will ever become popular (perhaps with the incoming Semantic web) or is it destined for academic work and system proofs only?
I don't believe Prolog will become any more popular than it is. The Japanese had some sort of a government project a few years ago, but I think that fell through. Now we have Erlang, Haskell and Java libraries that replicate most functionality, and then some.
It seems like I had to as part of an interview once. Still, I have needed to write recursive functions many, many times and it's good to review some of what can go wrong.
Thanks - I was struggling to remember where I had heard of a practical application of the Fibonacci sequence in CS. The Fibonacci Heap was it - I'm pretty sure it was mentioned on my CS course and it must have been '87, same year it was published!
Yeah, I recalled the name as random sort, and was pleasantly surprised when the bogosort link directed to the same algorithm. I particularly enjoyed the "Quantum Bogosort" algorithm on the Wiki page.
My personal favorite has always been permutation sort, where you try all possible permutations of a sequence and check if it is sorted.
What's nice about it is that it is deterministic yet ridiculously slow.
It can also be really easily implemented in Prolog[1] where you simply define what a permutation and being sorted means. After that you just search for a sorted permutation.
Dunno. IIRC Every time I have seen the "bad" Fibonacci recursive algorithm has been followed by the "good" recursive one (bottom-up), which is O(1) in size if your language/implementation does tail-call elimination...
The key point is that he is trying to calculate arbitrarily large fibonacci numbers, and so the addition is itself an O(n) operation. This means that even after memoization, the complexity is still O(n^2), and he uses a number of tricks to reduce that.
Bogosort terminates with probability 1. Is it really reasonable to say it doesn’t guarantee termination? People often do say that about randomised algorithms, which confuses me a little. Probability 1 is as guaranteed as anything probabilistic can reasonably hope to be, isn’t it?
[Edited to add: thanks for the replies. I think I expressed myself poorly here. It’s not that I don’t understand the difference between “with probability 1” and “always”. What I mean is that I don’t understand why people sometimes make a big deal out of it, and say “that algorithm is not even guaranteed to terminate” as though that somehow means the algorithm is untrustworthy or useless. The trouble with the game “roll a die until you get 100 sixes in a row” is that it has an astronomical expected running time – not that it might _never_ end, but that it will almost certainly take a very long time.]
Probability 1 is as guaranteed as anything probabilistic can reasonably hope to be, isn’t it?
Unfortunately, it isn't. This is the reason why in mathematics, we say that "possibility 1" means that an event is "almost sure", which is different from a "sure" event.
For instance, you can play a game where you roll a dice over an over again. You win if you get 100 times a 6 in sequence. If you play this game without any time limit, your possibility to win is exactly 1, which means that you can be "almost sure" you'll win in the end. However, you can't be sure, because there is still the possibility that you don't get 100 times a 6 in an eternity.
Note that this only happens in infinite probability spaces. In finite spaces, "almost sure" and "sure" are equivalent.
Also note that the same holds for "probability 0" which means "almost impossible", not to be confused with "impossible".
Bogosort terminates with probability 1. Is it really reasonable to say it doesn’t guarantee termination?
Yes.
Particularly when done on real machines that use pseudo-random numbers that will repeat in finite time. For instance many use http://en.wikipedia.org/wiki/Mersenne_twister with a version that will repeat after 2^19937-1 steps. With a random array of 2500 elements, the odds are very good that bogosort will fall into a loop before it has successfully sorted the array.
I think that depends on how you cycle the numbers. If you use the new positions as part of the randomization process then the average cycle time becomes much larger than 2^19937-1.
AKA: (Ignoring the fact it takes more than one number to shuffle a deck). After the first cycle of 2^19937-1 you have ~1 in (2500!) that you repeat the initial position. After the second cycle it's ~2 in (2500!) etc. It may be possible that your odds of success are around (2^19937-2)/(2^19937-1).
Of course how you cycle the numbers matters. I assumed that you take the original deck and keep on coming up with random shuffles of it. If so, then your cycle time will be the cycle time of the algorithm times (worst case) the length of the array. However if you do the obvious trick of doing each random shuffle to the previously shuffled deck, then the cycle time for the permutations to repeat becomes much, much longer - long enough that you are very, very likely to hit sorted at some point. But it is possible that you won't. And if you hit the case where you won't, then you'll absolutely never terminate.
That's an abuse of notation. Quicksort is expected O(N log N), where the expectation is usually over a uniform distribution of all permutations.
The actual notation for proportional asymptotic dominance of expectation is too convoluted for use, so people just don't use one. But to even leave out the word "expected" is just plain wrong.
In a somewhat similar genre, I enjoyed this paper on computing primes, and some of the very-inefficient algorithms that have been used for such purposes (often mislabeled as the "sieve of Eratosthenes"): http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
If you're writing in JavaScript it's a nice candidate for self memoizing functions. Obviously this is only helpful if you're making numerous requests to the method, a single request still produces a series of recursive calls.
Very good demonstration of subsequent improvements of a naive algorithm. To me that was somewhat depreciated by the fact that you can actually calculate n-the Fibonacci number using Binet's closed form formula (http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_ex...). You will need arbitrary precision arithmetic starting with certain 'n' though, as IEEE 754 will not give you correct result.
Given the precision needed in the arithmetic, I've always wondered if the matrix exponentiating method wouldn't be faster. I've never been motivated enough to actually research it, test it, or think hard about it though...
I’d love to see an exact algorithm using the closed form formula. It’s obviously possible to do — though it seems fairly complicated — but would it really be faster? My instinct is that it would be slower: if anyone wants to prove me wrong, I’d be thrilled!
(The asymptotic complexity is surely the same, in any case.)
You can’t get an exact answer in general by using fixed-size floating-point numbers, as your pseudocode does. The number of bits of precision you need to get an exact answer is going to depend on the input value, presumably linearly.
Note that you can optimize this by removing the "(1-c)^n" part, because its absolute value is smaller than 0.5 for big n, so it can't influence the (rounded) result.
This trick is especially funny when performed using a hand calculator. For instance, just calculate the powers of the golden ratio and observe numbers like "xxxxx,0000yyy" "xxxxx,9999yyy", which become closer and closer to integers.
> The closed form complexity is not the same (not sure how it could be)
I disagree. Both variants are more similar than it may seem at a first glance. In essence, the question here is which of the following calculations is more efficient for big n:
a) the power function for arbitrary precision floating point
b) the power function for 2x2 matrices of arbitrary precision integers
Assuming that we can ignore the initial effort to calculate the golden ratio (can we?), Binet's formula is based on a), while the article's best algorith is an application of b).
I'm absolutely undecided which one works better, mostly because both kinds of power functions are working essentially the same way, with a complexity of O(log(n)).
Also, the precisions required for both cases seem to be quite similar.
It's straightforward (if messy) to calculate two rationals exactly, then approximate √5 to within certain error bounds.
I wrote some code which does this, just as an experiment... and it's unbearably slow. Just try it.
import Control.Arrow
import Data.Ratio
newtype Matrix2 a = Matrix2 (a, a, a, a) deriving (Show, Eq)
instance (Num a) => Num (Matrix2 a) where
Matrix2 (a, b, c, d) * Matrix2 (e, f, g, h) =
Matrix2 (a*e+b*g, a*f+b*h, c*e+d*g, c*f+d*h)
fromInteger x = let y = fromInteger x in Matrix2 (y, 0, 0, y)
fib n = let Matrix2 (_, x, _, _) = Matrix2 (1, 1, 1, 0) ^ n in x
binom n =
scanl (\a (b, c)-> a*b `div` c) 1 $
takeWhile ((/=) 0 . fst) $ iterate (pred *** succ) (n, 1)
evens (x:_:xs) = x : evens xs
evens xs = xs
odds (_:x:xs) = x : odds xs
odds _ = []
iterate' f x = x : (iterate f $! f x)
powers b = iterate (b *) 1
esqrt e n = x where
(_, x):_ = dropWhile ((<=) e . abs . uncurry (-)) $ zip trials (tail trials)
trials = iterate (\x -> (x + n / x) / 2) n
fib' n = (a, b) where
d = 2 ^ n
a = sum (zipWith (*) (odds $ binom n) (powers 5)) % d
b = sum (zipWith (*) (evens $ binom n) (powers 5)) % d
fib2 n = numerator r `div` denominator r where
(a, b) = fib' n
l = lcm (denominator a) (denominator a)
r = a + esqrt (1 % l) (b * b / 5) + 1 % 2
I have not profiled the code, but I expect that the Newtonian approximation of sqrt is the slowest part.
It would surprise me if "powers" was a big problem relative to the rest. And no shortcuts: each of [1, b b², b³, b⁴, b⁵, …] (up to n) is used in the binomial expansion.
Hmm, I screwed up and wrote iterate instead of iterate' where I intended. Not that it really makes a difference.
In any case, the average time to print the first 20 Fibonacci numbers is 0.4157ms with matrix exponentiation and 1.707063s with exact integer math (yes, the units are different) on a T400.
Actually, you need just Z ring enriched (if this is right word) with sqrt(5). So instead of one float, you use pair of integers of arbitrary precision. So
This is a great method! It avoids arbitrary floating point arithmetics[1] and thus makes it more comparable with the method presented in the article.
As far as I can see, this algorithm is very similar to the matrix operations of the article, at least in terms of complexity.
[1] Only the very last step might involve fp arithmetics, because you'll have to convert the result pair (x,y) back to x+y*sqrt(5). However, we already know that the result has to be an integer, so the its second part will always be 0, no fp arithmetics required.
Well, both have computed the Fibonacci terms at a given level, so how different could it be? Here's my implementation:
def fib_fast2(n):
assert n >= 0
a, b = 2, 0 # invariant: a,b are components of 2(phi^n)
for bit in bits(n):
a, b = (a*a + 5*b*b)>>1, a*b
if bit: a, b = (a + 5*b)>>1, (a+b)>>1
return b
It's almost identical runtime as the one in the article - a hair slower (15.32s vs. 16.17s to compute fib 10M). They're probably related by some well known relation between Fibonacci numbers.
That’s very nice! And a little rearrangement will bring it from three to two “big multiplications”, making it faster than my routine:
def fib_ring2(n):
assert n >= 0
a, b = 2, 0 # invariant: phi^n = (a + b*sqrt(5)) / 2
for bit in bits(n):
ab = a*b
a, b = (a+b)*((a+5*b)//2) - 3*ab, ab
if bit: a, b = (a + 5*b)//2, (a+b)//2
return b
That method is, for me at least, really really fast.
Eleven seconds to compute the 1,000,000th Fibonacci number. Compared with the naive method which takes 143 seconds.
Code below:
def add(a, b): return (a[0]+b[0], a[1]+b[1])
def sub(a, b): return (a[0]-b[0], a[1]-b[1])
def divsq5(a): return (a[1], a[0]/5)
def mul(a, b): return (a[0]*b[0]+5*a[1]*b[1],a[0]*b[1]+a[1]*b[0])
def subi(x, a): return (x-a[0],-a[1])
def pow(b, e):
r = (1, 0)
while e > 0:
if e & 1 == 1:
result = mul(r, b)
e = e >> 1
b = mul(b, b)
return r
twophi = (1,1)
def fib0(n):
return divsq5(sub(pow(twophi, n), pow(subi(2,twophi), n)))[0]>>n
def fib1(n):
a = 0
b = 1
while n > 0:
(a,b) = (b, a+b)
n -= 1
return a
1. Given how your big integer code is probably implemented, you probably want to distribute the divisions by 2 into the intermediate calculations, rather just shifting by n at the very end.
2. Consider that (a, -b)(c, -d) = (ac + bd, -(bc + ad)). That means that the powers of (1,1) and (1,-1) are always related by negating the second component. Thus you only need to actually compute the power of (1,1) and double it.
Yeah, there are a whole lot of things I could do, but I just wanted a really quick program to see what the relative speeds were. And I have to say I was impressed.
I did this once for a Project Euler question. I called the function "awful-thing". It worked.
;this calculates (a+b√sqr)^n, using exponentiation by squaring
(def awful-thing (n a b sqr)
(afnwith (n n i a r b it 1 rt 0)
(if (is n 0)
(list it rt)
(even n)
(self (/ n 2)
(+ (square i) (* sqr (square r)))
(* 2 i r)
it
rt)
(self (- n 1)
i
r
(+ (* i it) (* sqr r rt))
(+ (* i rt) (* r it))))))
(def fib (n)
(with (pow2 (expt 2 n)
hh (map - (awful-thing n 1 1 5) (awful-thing n 1 -1 5)))
;now hh = '(a b) where (fib n) = (a + b*root5)/(root5*pow2)
;a should = 0
(/ (cadr hh)
pow2)))
I think it’s only the jargon that’s confusing you. As long as you can remember (or look up) the definition of matrix multiplication, then it genuinely is easy to verify.
I always find it incredibly difficult to concentrate on comparisons of Fibonacci sequence algorithms when I know for a fact that there is a closed-form expression[1] which gives F(n) in constant time.
72 comments
[ 3.1 ms ] story [ 93.5 ms ] threadThe code below consists of declaring that S is a sorted version of L as long as S is a permutation of L and S is sorted.
Where is ZS defined? Bear in mind that I only took half a module's worth (6 modules in a year at my university) of Prolog and we really just used it for learning predicate logic.
Nevertheless, using Prolog (especially when my code worked!) blew my mind. In the assignment I had to write a program that mimicked part of an aircraft controller in that each aircraft in it's flightplan had to be separated from all other aircraft. Do you think Prolog will ever become popular (perhaps with the incoming Semantic web) or is it destined for academic work and system proofs only?
I don't believe Prolog will become any more popular than it is. The Japanese had some sort of a government project a few years ago, but I think that fell through. Now we have Erlang, Haskell and Java libraries that replicate most functionality, and then some.
[Edit:] But there are some interesting data structures based on them: http://en.wikipedia.org/wiki/Fibonacci_heap
http://en.wikipedia.org/wiki/Bogosort#Quantum_bogosort
http://en.wikipedia.org/wiki/Quarantine_%28Greg_Egan_novel%2...
What's nice about it is that it is deterministic yet ridiculously slow.
It can also be really easily implemented in Prolog[1] where you simply define what a permutation and being sorted means. After that you just search for a sorted permutation.
[1] http://rosettacode.org/wiki/Sorting_algorithms/Permutation_s...
http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html...
The key point is that he is trying to calculate arbitrarily large fibonacci numbers, and so the addition is itself an O(n) operation. This means that even after memoization, the complexity is still O(n^2), and he uses a number of tricks to reduce that.
But it's a nice discussion in any case.
Nonono, Bogosort is way worse than naive recursive fibonacci - the former doesn't even guarantee termination, recursive fibonacci still does.
If you want to calculate fibonacci numbers not as a misguided exercise in algorithms but actually efficiently, use an algebraic form: http://en.wikipedia.org/wiki/Fibonacci_number#Computation_by...
[Edited to add: thanks for the replies. I think I expressed myself poorly here. It’s not that I don’t understand the difference between “with probability 1” and “always”. What I mean is that I don’t understand why people sometimes make a big deal out of it, and say “that algorithm is not even guaranteed to terminate” as though that somehow means the algorithm is untrustworthy or useless. The trouble with the game “roll a die until you get 100 sixes in a row” is that it has an astronomical expected running time – not that it might _never_ end, but that it will almost certainly take a very long time.]
http://en.wikipedia.org/wiki/Almost_surely
Bogosort almost surely terminates.
Unfortunately, it isn't. This is the reason why in mathematics, we say that "possibility 1" means that an event is "almost sure", which is different from a "sure" event.
For instance, you can play a game where you roll a dice over an over again. You win if you get 100 times a 6 in sequence. If you play this game without any time limit, your possibility to win is exactly 1, which means that you can be "almost sure" you'll win in the end. However, you can't be sure, because there is still the possibility that you don't get 100 times a 6 in an eternity.
Note that this only happens in infinite probability spaces. In finite spaces, "almost sure" and "sure" are equivalent.
Also note that the same holds for "probability 0" which means "almost impossible", not to be confused with "impossible".
The probability that the game ends is 1, but it is not guaranteed to end. Yet no reasonable person will worry about this in practice.
Your example with 100 times 6 in a row is conceptually exactly the same, just the expected time until the game terminates is much, much larger.
Yes.
Particularly when done on real machines that use pseudo-random numbers that will repeat in finite time. For instance many use http://en.wikipedia.org/wiki/Mersenne_twister with a version that will repeat after 2^19937-1 steps. With a random array of 2500 elements, the odds are very good that bogosort will fall into a loop before it has successfully sorted the array.
AKA: (Ignoring the fact it takes more than one number to shuffle a deck). After the first cycle of 2^19937-1 you have ~1 in (2500!) that you repeat the initial position. After the second cycle it's ~2 in (2500!) etc. It may be possible that your odds of success are around (2^19937-2)/(2^19937-1).
And the worst case of bogosort is not to terminate.
The actual notation for proportional asymptotic dominance of expectation is too convoluted for use, so people just don't use one. But to even leave out the word "expected" is just plain wrong.
(The asymptotic complexity is surely the same, in any case.)
The closed form complexity is not the same (not sure how it could be). Here's the code written in standard most language psuedocode:
Note: I assume your point wasn't that computing the golden ratio exactly is... well ummm... time consuming.You can’t get an exact answer in general by using fixed-size floating-point numbers, as your pseudocode does. The number of bits of precision you need to get an exact answer is going to depend on the input value, presumably linearly.
The answer should be:
This trick is especially funny when performed using a hand calculator. For instance, just calculate the powers of the golden ratio and observe numbers like "xxxxx,0000yyy" "xxxxx,9999yyy", which become closer and closer to integers.
I disagree. Both variants are more similar than it may seem at a first glance. In essence, the question here is which of the following calculations is more efficient for big n:
a) the power function for arbitrary precision floating point
b) the power function for 2x2 matrices of arbitrary precision integers
Assuming that we can ignore the initial effort to calculate the golden ratio (can we?), Binet's formula is based on a), while the article's best algorith is an application of b).
I'm absolutely undecided which one works better, mostly because both kinds of power functions are working essentially the same way, with a complexity of O(log(n)).
Also, the precisions required for both cases seem to be quite similar.
It's straightforward (if messy) to calculate two rationals exactly, then approximate √5 to within certain error bounds.
I wrote some code which does this, just as an experiment... and it's unbearably slow. Just try it.
Also, the multiplication might be a bottleneck, there has been a lot of research in this area. (the original article mentions some)
It might be a good idea to use an arbitrary precision library like GMP, which provides super-fast multiplication and thus power() algorithms.
(It also has an implementation of Fib() which outperforms both approaches here, but that's not my point.)
It would surprise me if "powers" was a big problem relative to the rest. And no shortcuts: each of [1, b b², b³, b⁴, b⁵, …] (up to n) is used in the binomial expansion.
In any case, the average time to print the first 20 Fibonacci numbers is 0.4157ms with matrix exponentiation and 1.707063s with exact integer math (yes, the units are different) on a T400.
(a,b)+(c,d) = (a+c,b+d)
(a,b)/sqrt(5) = (b,a/5)
(a,b)(c,d) = (ac+5bd,ad+bc)
2phi = (1,1)
F(n) = (2phi^n - (2-2phi)^n)/(2^n*sqrt(5))
[edit] As ot said, the right word is "extended"
As far as I can see, this algorithm is very similar to the matrix operations of the article, at least in terms of complexity.
[1] Only the very last step might involve fp arithmetics, because you'll have to convert the result pair (x,y) back to x+y*sqrt(5). However, we already know that the result has to be an integer, so the its second part will always be 0, no fp arithmetics required.
I think "extended" is the right word (see http://en.wikipedia.org/wiki/Ring_extension)
It is commonly denoted as Z[sqrt(5)]
Although I agree that this whole stuff looks very similar, I'm not sure whether it really leads to the exact same calculations.
Eleven seconds to compute the 1,000,000th Fibonacci number. Compared with the naive method which takes 143 seconds.
Code below:
1. Given how your big integer code is probably implemented, you probably want to distribute the divisions by 2 into the intermediate calculations, rather just shifting by n at the very end.
2. Consider that (a, -b)(c, -d) = (ac + bd, -(bc + ad)). That means that the powers of (1,1) and (1,-1) are always related by negating the second component. Thus you only need to actually compute the power of (1,1) and double it.
Since the Fibonacci numbers are defined by a linear recurrence, we can express the recurrence as a matrix, and it’s easy to verify that...
It's been way too long since my math minor for me to understand that.
I think it’s only the jargon that’s confusing you. As long as you can remember (or look up) the definition of matrix multiplication, then it genuinely is easy to verify.
That's what I call a bad algorithm!
[1] http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_ex...