"The worst algorithm in the world?" | Hacker News Reader

"The worst algorithm in the world?"

bosker.wordpress.com

161 pointsmccutchen69 comments

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

(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"

> Z ring enriched (if this is right word) with sqrt(5)

I think "extended" is the right word (see http://en.wikipedia.org/wiki/Ring_extension)

It is commonly denoted as Z[sqrt(5)]

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.

And when you calculate this out, you internally wind up doing the same calculation as the matrix method.

How do you know that? Have you actually done the math?

Although I agree that this whole stuff looks very similar, I'm not sure whether it really leads to the exact same calculations.

If you compute the eigenvalue decomposition of

  A={{0,1},{1,1}}

i.e.

  A = V^-1 D V

where D is a diagonal matrix with the eigenvalues phi and 1/phi as the diagonals, then

  A^n = V^-1 D^n V

where D^n is diagonal with phi^n and 1/phi^n as the diagonals. This gives exactly the above well known formula.

Yes, I have actually done the math. On the surface it looks very different, but a lot of the same numbers show up in intermediate calculations.

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

Oops, I should have found that. I looked and thought there was some reason it wouldn't work. Interesting post, thanks for the fun.

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

Two optimizations you can apply:

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

Thank you, that's a nice observation

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

http://stackoverflow.com/questions/4327846/calculating-fibon...

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 think the biggest issue in this code in the naive implementation of "powers". It runs in O(b) while it should take at most O(log(b)).

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

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 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:

   fib(n) {
      double c = 1.6180339887;
      return round((power(c, n) - power(1-c, n))/sqrt(5));
   }

Note: I assume your point wasn't that computing the golden ratio exactly is... well ummm... time consuming.

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

Indeed. I thought the OP was referring to the recursive algorithm, not the final algorithm.

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.

Actually I think that roughly was my point. :-)

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.

When you round you'll always get the exact integer answer, even using the precision of the golden ratio I gave in the code.

That’s really not true. You don’t have to take my word for it: try it with fib(1000), say.

The answer should be:

  43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875

Oh, the old n=1000 trick. Well played.

can't wait to break this out next time it comes up in an interview!

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

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?

ZS is defined by select.

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’s not just bad in the way that Bubble sort is a bad sorting algorithm; it’s bad in the way that Bogosort is a bad sorting algorithm.

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

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".

Your example is misleading, because the same can be said about the game where you roll a dice over and over again, and you win if you get a 6.

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.

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.

There is a distinction between "surely" and "almost surely":

http://en.wikipedia.org/wiki/Almost_surely

Bogosort almost surely terminates.

If you talk about O(something(n)) run times there's an implicit agreement that you are talking about worst case behavior.

And the worst case of bogosort is not to terminate.

Not sure about that: qick sort is generally said to be O(NlogN), but its worst case is O(N^2) (the simple implementation anyway).

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.

[1] http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_ex...

Raising something to the power of n is not 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

It seems the author hasn't read SICP:

http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html...

I don't think you've read the article, actually.

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.

His way around that is to use the matrix recurrence relation, which is also there in SICP as an exercise.

But it's a nice discussion in any case.

Ppffftt, I wrote a recursive algorithm that calculates 2^n in O(2^n) time for an exam question.

That's what I call a bad algorithm!

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.

But the single request will still be faster. Try it.

I always liked the sort where you randomize the elements, check to see if they're sorted, and if not try again.

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.

[1] http://rosettacode.org/wiki/Sorting_algorithms/Permutation_s...

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.

http://en.wikipedia.org/wiki/Bogosort#Quantum_bogosort

If you enjoyed that you might like the novel Quarantine by Greg Egan:

http://en.wikipedia.org/wiki/Quarantine_%28Greg_Egan_novel%2...

I never understood why Bogosort made the check so efficient; it seems so silly. Fortunately, somebody solved this problem and made bogobogosort: http://www.dangermouse.net/esoteric/bogobogosort.html

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

This started out taking baby steps, then took a huge leap in complexity right here:

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.

Sorry.

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.

  [ fib(n-1) fib(n)   ] x [0 1]
  [ fib(n)   fib(n+1) ]   [1 1]
  
  = [ fib(n)   fib(n-1)+fib(n) ]
    [ fib(n+1) fib(n)+fib(n+1) ]
  
  = [ fib(n)   fib(n+1) ]
    [ fib(n+1) fib(n+2) ]

Interesting writeup, thanks. OT: For all its fame, has any of you ever needed to calculate Fibbonacci numbers in real life? I haven't.

And neither am I a member of the Fibonacci Association http://www.mathstat.dal.ca/fibonacci/ http://www.fq.math.ca/

[Edit:] But there are some interesting data structures based on them: http://en.wikipedia.org/wiki/Fibonacci_heap

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.

Reply on news.ycombinator.com