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.
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)))
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.)
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.
> 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.
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.
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 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.
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.
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.
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
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.
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.
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.
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.
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...
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.
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!
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.
