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]

Log In

or

Sign in with Github

Sign in with Twitter

or

Sign in with Github

Sign in with Twitter

Sign Up

or

Sign in with Github

Sign in with Twitter

Reset Password

Enter the email address associated with your account, and we'll email you a link to reset your password.

ab initio.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!