Mathematically, why does this SICP algorithm for the exponent of a number modulo another number work?

Question

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):

1. First, we call (expmod 5 3 3). Mathematically, this means that we're asking for (5^3) mod 3.
2. As the second parameter is odd, we compute (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?

Answer 1

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

Answer 2

(ab) mod n = [a (b mod n)] mod n

is also true.

Here's a proof by induction on a.

Base case: when a = 0, (0b) mod n = 0 mod n = [0 (b mod n)] mod n.

Inductive case:

By induction hypothesis, assume that (ab) mod n = [a (b mod n)] mod n is true. We need to prove that ((a+1) b) mod n = [(a + 1) (b mod n)] mod n.

((a+1) b) mod n
= (ab + b) mod n
= (ab mod n) + (b mod n)
= [a (b mod n)] mod n + (b mod n)             by induction hypothesis
= [a (b mod n)] mod n + (b mod n) mod n
= [a (b mod n) + (b mod n)] mod n
= [(a + 1) (b mod n)] mod n

as desired.

This concludes the proof that

(ab) mod n = [a (b mod n)] mod n

---

In fact, you can see that

(ab) mod n = [(a mod n) (b mod n)] mod n

is a result that follows from it. Here's a proof:

(ab) mod n 
= [a (b mod n)] mod n            by what we just proved
= [(b mod n) a] mod n
= [(b mod n) (a mod n)] mod n    by what we just proved
= [(a mod n) (b mod n)] mod n