Hacker Newsnew | past | comments | ask | show | jobs | submit | mschoenert's commentslogin

This prime generating code has fascinated me for a while.

It was first described by E.W. Dijkstra "Notes on Structured Programming (EWD249), 2nd; 1970; TH Eindhoven". https://www.cs.utexas.edu/~EWD/ewd02xx/EWD249.PDF

And then reformulated by D.E. Knuth "Literate Programming; 1984" http://www.literateprogramming.com/knuthweb.pdf

Both use it to demonstrate their way of deriving and documenting an algorithm.

And then R. Martin used it again in his book "Clean Code; 2008". (Though I'm not 100% certain what he wants to demonstrate with his rather difficult formulation.)

For the sake of this email I am going to call it "Dijkstra's algorithm".

If we look at it from an engineering perspective, then Dijkstra's algorithm is never the best choice (at least not in 2025).

If performance is not the most important aspect, e.g. if we need less than 100000 primes, then the straightforward trial division is just as fast - and soooo much simpler (basically impossible to get wrong). Note that this is different in 2025 than it was in 1970, back then multiplications and divisions were much more expensive so it made sense to minimize them.

If performance is the most important aspect, then the sieve of Eratosthenes is much faster. And you can implement the sieve with a sliding window, which is just as fast (or even a bit faster because of caching) and uses a bounded amount of memory.

Concretely - my implementation of the sliding window sieve of Eratosthenes is about 50 times faster than Dijkstra's algorithm (on the same machine) when generating the first 200 million primes (7.5 sec vs 6 min).

The reformulation - both by Ousterhout and Martin - computes multiples for every new prime found. This will quickly lead to overflow, e.g. the 6543-th prime is 65537, so its square will overflow 32 bit unsigned integers. Dijkstra's formulation on the other hand only computes multiples that are actually used to filter, the multiples are never (much) larger than the candidate.

Note that computing the multiples so eagerly will also make the multiples vector unnecessarily large (wasting memory).

Knuth and Dijkstra both remark that there is a subtle problem with the correctness of the algorithm. When incrementing lastMultiple and then later accessing primes[lastMultiple] und multiples[lastMultiple] it is not obvious that those will already be assigned. In fact they will be, but it is very difficult to prove that. It is a consequence of the fact that for every integer n there is always a prime between n and n*n, which follows from Betrand's theorem, a difficult number-theoretical result.

So if you look at Ousterhout's and Martin's reformulation and think "a yes - now I get the algorithm", then beware: You've missed at least two aspects that are relevant for the algorithm's correctness. ;-)


As pointed out above - the most efficient to compute large fibonacci numbers is to compute the matrix power [[1,0],[1,1]]^n using repeated squaring. Or you could use the known identities to compute Lucas numbers, which amount to the same thing. The lispmath talks about computing fib(40000) in 100..200 ms, the repeated squaring approach computes the same number in < 5 ms (on my not very powerful machine).


My thoughts at this point in the original article was:

Well - if you want to demonstrate how mathematics helps here, then you should mention that computing the n-th fibonnaci number can be computed as the n-th power of the matrix [[1,0],[1,1]], and that the most efficient way to compute powers (in any associative domain) is the repeated squaring algorithm.


There are many ways to obtain the formula mentioned in the article. One of them is by diagonalizing this matrix and then applying exponentiation. How can another algorithm be better than constant time in the most general case?


No - there isn't really an O(1) solution. And the reason is that the Fibonacci numbers grow without limit. So there can't be an O(1) algorithm - even just writing down the answer takes O(N) - because the answer has O(N) bits.

Concretely Fibonacci(1480) is the largest Fibonacci number that fits into a double. So for higher Fibonacci numbers you need to compute this with arbitrary size integers (or floats). And then look at it in terms of bit complexity (i.e. the time to compute a product grows with the size of the integers).

Put differently. The "constant time" algorithm only works for N up to 1480. And if you limit N, then it doesn't really make sense to talk about the asymptotic complexity.

I believe the optimal algorithm is to compute the Fibonacci number by computing the N-th power of the matrix [[1,1],[1,0]], and compute the power using the repeated squaring algorithm. There are alternative formulations of that using identities for Fibonacci or Lucas numbers, but those identities basically correspond to squaring the matrix.

If you really want the asymptotically fastest algorithm you must use the fastest algorithm for integer multiplication (Harvey and van der Hoeven’s).

But even with the naive multiplication for arbitrary size integers you can compute Fibonacci(1000000) - which is a number with more than 200'000 digits - in less than a second. ;-)


Well, with Scheme, the difference on computing fib(n) recursively vs the iterative way shows up really fast.


The idea that the "Lisp programmers" are somehow adverse to mathematics is historically untenable. Macsyma was written in Lisp. And Macsyma was one of - if not THE most important application for the Symbolic Lisp machines (the ancestor of all later Lisp machines). Many later systems (Maple, Mathematica, ...) were written by people who would consider themselves belonging to the Lisp crowd. For example all of those systems had automatic memory management with garbage collection, most had closures/lambdas, etc.


Guidelines | FAQ | Lists | API | Security | Legal | Apply to YC | Contact

Search: