schemesicpexponentiationmodular-arithmetic

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


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?


Solution

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