fast computation of Fibonacci numbers with quick matrix exponentiation (as we were taught to do it in school)

fib.rb

# 2x2 matrix multiplication
#  [a  b]       [e  f]     [a*e+b*g a*f+b*h]
#  [c  d]   *   [g  h]  =  [c*e+d*g c*f+d*h]
def mult m1, m2
  a, b, c, d = m1
  e, f, g, h = m2
  [a*e+b*g, a*f+b*h, c*e+d*g, c*f+d*h]
end

# quick exponentiation
# x^(2k)   = (x^k)^2
# x^(2k+1) = x*(x^k)^2 
def exp m, n
  if n == 1
    m
  elsif n % 2 == 0
    mult(o = exp(m, n/2), o)
  else
    mult(m, mult(o = exp(m, n/2), o))
  end
end


# fibonacci
#                           [0 1]
# [fn-1 fn] = [fn-2 fn-1] * [1 1]
#
#                       [0 1] 
#           = [f0 f1] * [1 1] ^ n
def fib n
  if n < 2
    [0, 1][n]
  else
    exp([0, 1, 1, 1], n).last
  end
end


# demo
(0..20).each{|i| puts "fib[#{i}] = #{fib(i)}"}