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: Program Derivation
