Fast Exponentiation

Today Shin and I were grading the final exam of Program Construction and Reasoning in FLOLAC 2008 and one of the students presented an interesting derivation for the O(lg n) fast exponentiation algorithm where n is the exponent. In fact, he derived two different solutions by induction and fold-fusion respectively, the second of which is the standard version. Shin later observed that the first solution (the interesting one) is a higher-order fold. Below I'll describe the solution in (more or less) point-free style, but of course it would be much easier to do the derivation by induction.

I prefer to define unfoldr on Maybe functors.

unfoldr :: (b -> Maybe (a, b)) -> b -> [a]
unfoldr f s = case f s of Nothing -> []
                          Just (x, s') -> x : unfoldr f s'

A function binary can be defined to produce a "reversed binary bit list" representing the input number.

binary :: Integer -> [Bool]
binary = unfoldr (\n -> if n == 0 then Nothing
                                  else Just (odd n, n `div` 2))

For example, binary 11 = [True, True, False, True]. Its left inverse decimal can be defined as a fold.

decimal :: [Bool] -> Integer
decimal = foldr (\x y -> if x then 1 + 2 * y else 2 * y) 0

Let exp m n = mn. A naïve implementation of exp requires n - 1 multiplications. But we can introduce a hylomorphism and do fold-fusion to reduce the number of multiplications to O(lg n):

  exp m n
=   { identity }
  (exp m . id) n
=   { decimal is a left inverse of binary }
  (exp m . decimal . binary) n
=   { fusion, see below }
  (binexp m . binary) n

where binexp m . binary takes only O(lg n) time, since the length of the intermediate list is only O(lg n) and both binexp m and binary runs in time linear to the length of the intermediate list. The function binexp is specified by

binexp m xs = exp m (decimal xs)

The standard solution can be derived by directly fusing exp m into decimal, which gives

binexp m []           = 1
binexp m (True  : xs) = m * square (binexp m xs)
binexp m (False : xs) = square (binexp m xs)

The interesting solution, on the other hand, may be derived by first introducing flip on both sides:

flip binexp xs m = flip exp (decimal xs) m

After some eta-conversions we would get a point-free specification of flip binexp:

flip binexp = flip exp . decimal

Now it's obvious that we should try to fuse flip exp into decimal. The fusion condition is justified as follows:

Therefore flip binexp = foldr g (flip exp 0) = foldr g (const 1). (In fact this involves the debate of whether 0⁰ is defined. Currently I am convinced by Knuth's argument --- see Two Notes on Notation.) Applying flip to both sides again, we get binexp = flip (foldr g (const 1)), which expands to

binexp m []           = 1
binexp m (True  : xs) = m * binexp (square m) xs
binexp m (False : xs) = binexp (square m) xs

Notice that the last equation is a tail-recursive one.

--
Derivations using higher-order fusion always look much more frightening than the corresponding pointwise, inductive derivations...