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.