Because the exponentiation method lets you compute an arbitrary element of the Fibonacci sequence without calculating most of the ones before it. Your `benchmark' calculates them in ascending order, so obviously a straightforward memoizing/iterative implementation does well. Instead, it should just calculate the millionth Fibonacci number (say), ab initio.

You're doing too much work in the matrix multiplication: the matrices are always symmetrical, but you're not taking advantage of that. But there's a slicker way, which involves some slightly fancier mathematics.

Suppose that we know that t^2 = t + 1 for some t. (Equivalently, we'll work in the quotient ring R = Z[t]/(t^2 — t — 1).) Now, let's look at the powers of t. By using the rule that t^2 = t + 1, we will end up with only units and coefficients of t. (That this is the characteristic polynomial of the matrix discussed in the main article is, of course, not a coincidence.)

t^0 = 1
t^1 = t
t^2 = t + 1
t^3 = t^2 + t = (t + 1) + t = 2 t + 1
t^4 = 2 t^2 + t = 2 (t + 1) + t = 3 t + 2
t^5 = 3 t^2 + 2 t = 3 (t + 1) + 2 t = 5 t + 3

Hmm. There's a pattern: t (a t + b) = a t^2 + b t = a (t + 1) + b t = (a + b) t + a. I feel a theorem coming on: t^k = F_k t + F_{k-1} for all integers k. Proof. (a) t^1 = t = F_1 t + F_0. (b) Suppose the claim holds for k; then t^{k+1} = t t^k = (F_k + F_{k-1}) t + F_k = F_{k+1} t + F_k. (c) Notice that t (t — 1) = t^2 — t = (t + 1) — t = 1, so t^{-1} = t — 1; then, supposing again that the claim holds for k, we have t^{k-1} = (t — 1) (F_k t + F_{k-1}) = F^k t^2 + F_{k-1} t — F_k t — F_{k-1} = F_k (t + 1) + F_{k-1} t — F_k t — F_{k-1} = F_{k-1} t + (F_k — F_{k-1}) = F_{k-1} t + F_{k-2}, as required. []

Now, we can compute powers of t using the ordinary square-and-multiply method (or fancier methods involving windowing, combs, or whatever). The important pieces are:

(a t + b)^2 = a^2 (t + 1) + 2 a b t + b^2 = (2 a b + a^2) t + (a^2 + b^2) t
(a t + b) (c t + d) = a c (t + 1) + (a d + b c) t + b d = (a d + a c + b c) t + (a c + b d)

Here's a simple left-to-right exponentiation algorithm. Suppose we want to compute x^n, and we know that n = 2^k u + v, for 0 <= v < 2^k, and x^u = y. (Initially, we can arrange this with 2^k > n, so u = 0, v = n, and x = 1.) Now, suppose that v = 2^{k-1} b + v', where b is either zero or one, and v' < 2^{k-1}. Let u' = 2 u + b. Then n = 2^{k-1} u' + v', and x^{u'} = x^{2u+b} = y^2 x^b. Continue until k = 0.

Annoyingly, this means we need to count the number of bits in n, i.e., k = L(n), with 2^{i-1} <= n < 2^i. We can start with an upper bound, finding the smallest 2^{2^j} > n. Suppose, inductively, that we know how to determine L(m) for 0 <= m < 2^i, and that n < 2^{2i}. If n < 2^i, then we're done; otherwise write n = 2^i n' + r where 0 <= r < 2^i. Since n < 2^{2i}, we know that n' <= n/2^i < 2^i, so we can determine k' = L(n'). Let k = k' + i: I claim that 2^{k-1} <= n = 2^i n' + r < 2^k. Taking the lower bound first, we're given that 2^{k'-1} <= n', so 2^{k-1} = 2^i 2^{k'-1} <= 2^i n' <= 2^i n' + r. But also, n' < 2^{k'}, so n' + 1 <= 2^{k'} and 2^i n' + 2^i <= 2^i 2^{k'} = 2^k. But r < 2^i, so n = 2^i n' + r < 2^k as required. []

Enough theory! Now let's turn this into actual code!

#! /usr/bin/python

def power(x, n, id, sqr, mul):
  """
  Compute and return X^N for nonnegative N.

  This is an abstract operation, assuming only that X is an element of a
  monoid: ID is the multiplicative identity, and the necessary operations are
  provided as functions: SQR(Y) returns Y^2, and MUL(Y, Z) returns the
  product Y Z.
  """

  ## First, we need to find a 2^{2^i} > n.
  i = 1
  while (n >> i) > 0: i = 2*i

  ## Refine our estimate, so that 2^{k-1} <= n < 2^k.
  k, m = 1, n
  while m > 1:
    i = i >> 1
    mm = m >> i
    if mm > 0: m, k = mm, k + i

  ## Now do the square-and-multiply thing.
  y = id
  for i in xrange(k - 1, -1, -1):
    y = sqr(y)
    if (n >> i)%2: y = mul(x, y)

  ## We're done.
  return y

def fibonacci(n):
  """
  Compute and return the Nth Fibonacci number.

  Let F_0 = 0, F_1 = 1, F_{k+1} = F_k + F_{k-1}, for all integers k.  Note
  that this is well-defined even for negative k.  We return F_N.
  """

  ## Auxiliary functions for working in our polynomial ring.
  def poly_sqr((a, b)):
    a2 = a*a
    return 2*a*b + a2, a2 + b*b
  def poly_mul((a, b), (c, d)):
    ac = a*c
    return a*d + b*c + ac, ac + b*d

  ## Do the job.  For negative indices, we take powers of t^{-1}.
  if n < 0: return power((1, -1), -n, (0, 1), poly_sqr, poly_mul)
  else: return power((1, 0), n, (0, 1), poly_sqr, poly_mul)

## Command-line interface.
if __name__ == '__main__':
  import sys
  try: n = (lambda _, n: int(n))(*sys.argv)
  except TypeError, ValueError:
    print >>sys.stderr, 'Usage: fib N'
    sys.exit(1)
  print fibonacci(n)[0]