### 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) 0Let

`exp m n = m`^{n}

. 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) nwhere

`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) mAfter some eta-conversions we would get a point-free specification of

`flip binexp`

:
flip binexp = flip exp . decimalNow 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^{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) xsNotice 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

<< 回到主頁