$\newcommand{\defeq}{\mathrel{\mathop:}=}$

## 2008/07/11

### 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 00 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...

Labels: