/*****************************************************************************
 * File: ZeckendorfLogarithm.hh
 * Author: Keith Schwarz (htiek@cs.stanford.edu)
 *
 * An implementation of a function for computing logarithms based on the
 * Zeckendorf decomposition.  Zeckendorf's theorem states that any number can
 * be written as the sum of unique Fibonacci numbers such that no two
 * consecutive Fibonacci numbers are used.  For example:
 *
 *   15 = 13 + 2
 *   95 = 65 + 21 + 8 + 1
 *
 * The proof of Zeckendorf's theorem is actually quite simple and also gives
 * an algorithm for computing such a decomposition.  The proof is by
 * induction.  As base cases, both 0 and 1 are Fibonacci numbers.  For the
 * inductive step, assume that for all n' < n the claim holds and consider n.
 * If n is a Fibonacci number then the claim is trivially true.  Otherwise,
 * let F(k) be the largest Fibonacci number smaller than n and consider the
 * number n - F(k).  By the inductive hypothesis, n - F(k) can be written as
 * the sum of unique Fibonacci numbers with no two consecutive numbers used.
 * Consider the sum of those unique numbers, plus F(k).  If we can show that
 * F(k) doesn't appear twice, we have that each number is unique, and if we
 * can show that F(k-1) doesn't appear anywhere in the sum, then we have that
 * no two consecutive Fibonacci numbers are used.
 *
 * To see that F(k) doesn't appear twice, assume for the sake of contradiction
 * that it does.  This means that F(k) must be in the decomposition of 
 * n - F(k), and so we have that n - F(k) <= F(k), or that n <= 2F(k).  But
 * F(k) = F(k-1) + F(k-2) (assuming that F(k) >= 2; if it's not, then n = 0 or
 * n = 1 and we already know that the claim holds).  This means that
 * n <= F(k) + F(k-1) + F(k-2) = F(k+1) + F(k-2), contradicting the fact that
 * F(k) is the largest Fibonacci number smaller than n.
 *
 * To see that F(k-1) doesn't appear in the decomposition of n - F(k) we use
 * a similar proof by contradiction.  Assume that F(k-1) does appear in the
 * sum; then we have that n - F(k) <= F(k-1), so n <= F(k) + F(k-1) = F(k+1),
 * again contradicting that F(k) is the largest Fibonacci number less than n.
 *
 * Notice that this proof gives us a very simple algorithm for computing the
 * Zeckendorf decomposition of a number - just continuously subtract out the
 * largest Fibonacci number less than n and add it to the result.  It is
 * possible to compute this in O(lg n) time and O(1) space using the following
 * observation:
 *
 *   Given two consecutive Fibonacci numbers F(k) and F(k+1), we can compute
 *   the pairs (F(k-1), F(k)) and (F(k+1), F(k+2)) in constant time.
 *
 * This is true because
 *
 *   F(k-1) = F(k+1) - F(k)
 *
 * and
 *
 *   F(k+2) = F(k+1) + F(k)
 *
 * In other words, we can think of a pair of adjacent Fibonacci numbers as a
 * sort of "bidirectional iterator" that can advance forwards and backwards
 * one step in the Fibonacci sequence.
 *
 * Given such an iterator, we can compute the Zeckendorf decomposition in
 * O(lg n) time, O(1) space using the following trick:
 *
 *   1. Begin a "Fibonacci iterator" at (0, 1).
 *   2. Until the second term of the iterator exceeds n, advance the iterator.
 *   3. While n is not zero:
 *      1. If the first term of the iterator is less than n, add the first 
 *         term of the iterator to the resulting decomposition, then subtract
 *         the first term from n.
 *      2. Back the iterator up a step.
 *
 * This uses only O(1) memory (plus any memory used to store the resulting
 * sequence).  To see that it runs in O(lg n) time, recall that
 *
 *    F(k) = Theta(Phi^k)
 *
 * where Phi = (1 + sqrt(5)) / 2 is the Golden Ratio.
 *
 * The power of using this decomposition of a number is that the Fibonacci
 * sequence can be generalized to arbitrary groups, not just the natural
 * numbers under addition.  In fact, given any group G with generator g, we
 * can define the generalized Fibonacci sequence on G as
 *
 *    F_G (0)   = g^0 = 1
 *    F_G (1)   = g^1 = g
 *    F_G (n+2) = F_G(n)F_G(n+1) = g^(F(n) + F(n+1))
 *
 * (That last equality can be verified by a quick inductive argument)
 *
 * This provides us an extremely elegant algorithm for computing logarithms
 * using only O(lg lg n) arithmetic operations and O(1) space.  The idea is as
 * follows.  Suppose that we are given n = a^b and want to compute log_a(n).
 * Then we can compute b by computing its Zeckendorf decomposition in the
 * generalized Fibonacci with generator a.  More specifically:
 *
 *   1. Begin a "generalized Fibonacci iterator" at (a^0, a^1).  Record the
 *      exponents of both terms in the iterator at all times.
 *   2. Until the second term of the iterator exceeds n, advance the iterator.
 *   3. While n is not zero:
 *      1. If the first term of the iterator is less than n, add the
 *         exponent of the first term to the result, then divide n by the
 *         first term in the iterator.
 *      2. Back the iterator up a step.
 *
 * Assuming that b is integral, then this solution will yield an exact answer
 * after doing O(lg b) arithmetic operations.  Since b = O(lg n), this
 * algorithm runs using only O(lg lg n) arithmetic operations!  Moreover, it
 * uses memory for only a constant number of integers.
 */

         
#ifndef ZeckendorfLogarithm_Included
#define ZeckendorfLogarithm_Included

#include <utility>   // For pair
#include <stdexcept> // For std::domain_error

/**
 * Function: Integer ZeckendorfLogarithm(T n, const T& a);
 * Usage: cout << ZeckendorfLogarithm(256, 2) << endl; // Prints 8
 * ---------------------------------------------------------------------------
 * Given two values n and a, returns log_a n.  It is assumed that n is
 * an integral power of a; if not, the function rounds down.  If the base or
 * exponent is zero, a domain_error exception is thrown.
 *
 * The parameters must be of a type that supports multiplication and division.
 * It also must allow for the numeric values 0 and 1 to be cast to an element
 * of its type.  The elements of type T must be totally ordered and 
 * comparable.
 */

template <typename Integer, typename T>
Integer ZeckendorfLogarithm(T n, const T& a) {
  /* Check that the input is valid. */
  if (n == T(0))
    throw std::domain_error("Cannot take the log of zero.");
  if (a == T(0) || a == T(1))
    throw std::domain_error("Cannot take log base zero or one.");

  /* Construct our "Fibonacci iterator" to keep track of where we are in the
   * Fibonacci sequence, including both the base and the logarithm.
   */

  std::pair<T, T>             fibExponentIterator(T(1), a);
  std::pair<Integer, Integer> fibBaseIterator(0u1u);

  /* Keep marching the iterator forward until we arrive at n. */
  while (fibExponentIterator.second < n) {
    /* Compute the next value in the generalized Fibonacci sequence as
     *
     *   (F(k), F(k+1)) -> (F(k+1), F(k+2))
     */

    T nextExp = fibExponentIterator.first * fibExponentIterator.second;
    fibExponentIterator.first = fibExponentIterator.second;
    fibExponentIterator.second = nextExp;

    Integer nextBase = fibBaseIterator.first + fibBaseIterator.second;
    fibBaseIterator.first = fibBaseIterator.second;
    fibBaseIterator.second = nextBase;
  }

  Integer result = Integer(0);

  /* Keep dividing out terms from the Zeckendorf representation of the number
   * until it reaches one.
   */

  while (fibBaseIterator.first != Integer(0)) {
    /* If the first of the Fibonacci numbers is less than the current number,
     * divide it out and add the corresponding Fibonacci number to the
     * total.
     */

    if (fibExponentIterator.first <= n) {
      n /= fibExponentIterator.first;
      result += fibBaseIterator.first;
    }

    /* Walk both iterators backward a step:
     *
     *   (F(k), F(k+1)) -> (F(k-1), F(k))
     */

    T prevExp = fibExponentIterator.second / fibExponentIterator.first;
    fibExponentIterator.second = fibExponentIterator.first;
    fibExponentIterator.first  = prevExp;

    Integer prevBase = fibBaseIterator.second - fibBaseIterator.first;
    fibBaseIterator.second = fibBaseIterator.first;
    fibBaseIterator.first = prevBase;
  }

  return result;
}

#endif