O(log n) fibonacci numbers

fibonacci.cpp

#include <utility>

// The 2x2 fibonacci matrix
//
// { F(n+1) F(n)   }
// { F(n)   F(n-1) }
//
// can be represented as just the bottom row, since the top row can be computed:
// so just use a pair.
template <typename N>
struct multiply_fib
{
  constexpr std::pair<N, N> operator()(const std::pair<N, N>& x,
                                       const std::pair<N, N>& y) noexcept
  {
    return { x.first * (y.first + y.second) + x.second * y.first,
             x.first * y.first + x.second * y.second };
  }
};


// (Bottom row of) 2x2 matrix identity { 1, 0, 0, 1 }
template <typename N>
constexpr std::pair<N, N> identity_element(const multiply_fib<N>&) noexcept
{
  return { N(0), N(1) };
}

// Generic power function (see Elements of Programming or From Mathematics to
// Generic Programming)
template <typename A, typename N, typename Op>
constexpr A power(A a, N n, Op op) noexcept
{
  if (n == 0) return identity_element(op);

  while ((n & 1) == 0)
  {
    a = op(a, a);
    n >>= 1;
  }

  A r = a;
  n >>= 1;
  while (n != 0)
  {
    a = op(a, a);
    if ((n & 1) != 0)
    {
      r = op(r, a);
    }
    n >>= 1;
  }
  return r;
}

// The nth fibonacci number is the bottom left of the matrix { 1, 1, 1, 0 }
// raised to the power n
template <typename R, typename N>
R fibonacci(N n)
{
  if (n == 0) return R(0);
  return power(std::pair{ R{1}, R{0} },
               n,
               multiply_fib<R>()).first;
}

/*
#include <iostream>

int main()
{
  // produces 70th fibonacci number: 190392490709135
  std::cout << fibonacci<uint64_t, uint64_t>(70) << std::endl;
  return 0;
}
*/