Section 1.2.6 of SICP gives the following procedure:
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (square (expmod base (/ exp 2) m))
m))
(else
(remainder (* base (expmod base (- exp 1) m))
m))))
The authors claim that it "computes the exponential of a number modulo another number". For example (expmod 5 3 n)
should return (5^3) mod n.
However, from a mathematical point of view, I just can't see how it works. As reinforced by footnote 46, it is intended to use the property that for any positive integers a, b, and n, (ab) mod n = [(a mod n)(b mod n)] mod n, but I fail to see how it is actually using it. Consider (expmod 5 3 3)
:
(expmod 5 3 3)
. Mathematically, this means that we're asking for (5^3) mod 3.(remainder (* 5 (expmod 5 (- 3 1) 3)) 3)
i.e. (remainder (* 5 (expmod 5 2 3)) 3)
. Mathematically, this is [5 * [(5^2) mod 3]] mod 3. As the initial 5 in this expression does not have a mod 3 attached, this expression is not in the (ab) mod n = [(a mod n)(b mod n)] mod n form, so it fails to use the intended property.So, given that this it appears to not be using the intended property, why does this algorithm work? What property of modular arithmetic have I overlooked?
This is the definition of fast-exp from 1.2.4 that is refered to:
(define (fast-expt b n)
(cond ((= n 0) 1)
((even? n) (square (fast-expt b (/ n 2))))
(else (* b (fast-expt b (- n 1))))))
If we rename things to closer match expmod it looks like this:
(define (expt base exp)
(cond ((= exp 0) 1)
((even? exp)
(square (expt base (/ exp 2))))
(else
(* base (expt base (- exp 1))))))
To get a naive expmod
we can, for now, just calculate the remainder of each clause:
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (square (expt base (/ exp 2))) m))
(else
(remainder (* base (expt base (- exp 1))) m))
So far we have not used the footnote (ab) mod m = ((a mod m)(b mod m) mod m)
. Of course a special case of this is (aa) mod m = ((a mod m)(a mod m) mod m)
which gives (remainder (square a) m) = (remainder (sqaure (remainder a m)) m)
. We can use this with the even
clause so that
(remainder (square (expt base (/ exp 2))) m)
becomes:
(remainder (square (remainder (expt base (/ exp 2)) m))
m)
In the middle of this we have the remainder of an exponent so this is equivalent to:
(remainder (square (expmod base (/ exp 2) m)) m)
Using the new even
clause we have
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (square (expmod base (/ exp 2) m))
m))
(else
(remainder (* base (expt base (- exp 1))) m))
To simplify the odd clause lets use E
in place of (expt base (- exp 1))
for now.
By using the defining properties of mod
we can say for any number a
:
a = (+ (* (quotient a m) m) (remainder a m))
So it's also true that:
E = (+ (* (quotient E m) m) (remainder E m))
substituting this into our odd
clause:
(remainder (* base E) m)
gives:
(remainder (* base (+ (* (quotient E m) m) (remainder E m))) m)
We can ignore (* (quotient E m) m)
because any term containg this is divisible by m
and so will evaluate to 0
when doing the outer remainder
, so this equivalent to:
(remainder (* base (remainder E m)) m)
expanding E to it's orignal value:
(remainder (* base (remainder (expt base (- exp 1)) m)) m)
once again, in the middle, we have the remainder of an exponent so this becomes:
(remainder (* base (expmod base (- exp 1) m)) m)
And our expmod
is now:
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (square (expmod base (/ exp 2) m))
m))
(else
(remainder (* base (expmod base (- exp 1) m))
m))))