title: SICP Exercise 1.28 - Miller-Rabin - "at least half the numbers will reveal a nontrivial square root of 1 modulo n"
body: <p><strong>SICP Exercise 1.28</strong></p>
  
  <p><a href="https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#%_thm_1.28" rel="nofollow noreferrer">https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-11.html#%_thm_1.28</a></p>
  
  <blockquote>
    <p>One variant of the Fermat test that cannot be fooled
    is called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts
    from an alternate form of Fermat's Little Theorem, which states that
    if n is a prime number and a is any positive integer less than n, then
    a raised to the (n - 1)st power is congruent to 1 modulo n. To test
    the primality of a number n by the Miller-Rabin test, we pick a random
    number a &lt; n and raise a to the (n - 1)st power modulo n using the
    expmod procedure. However, whenever we perform the squaring step in
    expmod, we check to see if we have discovered a ``nontrivial square
    root of 1 modulo n,'' that is, a number not equal to 1 or n - 1 whose
    square is equal to 1 modulo n. It is possible to prove that if such a
    nontrivial square root of 1 exists, then n is not prime. It is also
    possible to prove that <strong>if n is an odd number that is not prime, then,
    for at least half the numbers a &lt; n, computing a^(n-1) in this way will
    reveal a nontrivial square root of 1 modulo n.</strong> (This is why the
    Miller-Rabin test cannot be fooled.) Modify the expmod procedure to
    signal if it discovers a nontrivial square root of 1, and use this to
    implement the Miller-Rabin test with a procedure analogous to
    fermat-test. Check your procedure by testing various known primes and
    non-primes. Hint: One convenient way to make expmod signal is to have
    it return 0.</p>
  </blockquote>
  
  <p>I've written my own solution and its results are consistent with the solutions provided here:</p>
  
  <p><a href="http://community.schemewiki.org/?sicp-ex-1.28" rel="nofollow noreferrer">http://community.schemewiki.org/?sicp-ex-1.28</a></p>
  
  <p><code>15</code> is an odd number that is not prime, so for at least half the numbers <code>a</code> from <code>1</code> to <code>14</code>, I expect computing <code>expmod(a, 14, 15)</code> will reveal a nontrivial square root of 1 modulo n, which is signified by <code>expmod</code> returning <code>0</code>.</p>
  
  <p>However, these are the results I get:</p>
  
  <pre><code>(expmod 1 14 15)
  &gt; 1
  (expmod 2 14 15)
  &gt; 4
  (expmod 3 14 15)
  &gt; 9
  (expmod 4 14 15)
  &gt; 0
  (expmod 5 14 15)
  &gt; 10
  (expmod 6 14 15)
  &gt; 6
  (expmod 7 14 15)
  &gt; 4
  (expmod 8 14 15)
  &gt; 4
  (expmod 9 14 15)
  &gt; 6
  (expmod 10 14 15)
  &gt; 10
  (expmod 11 14 15)
  &gt; 0
  (expmod 12 14 15)
  &gt; 9
  (expmod 13 14 15)
  &gt; 4
  (expmod 14 14 15)
  &gt; 1
  </code></pre>
  
  <p>As can be seen, only 2 of these results are <code>0</code>, which is way short of at least 7 as expected.</p>
  
  <p>Am I misunderstanding the statement? Am I being a complete idiot? Is the code wrong? Is SICP wrong? Many thanks.</p>
  
  <p>Edit 1: it was requested that I supply the exact code I'm using. Here it is, although I'm essentially just copying the solution I linked to, and aliasing <code>remainder</code> as <code>mod</code> because that's what my interpreter calls it.</p>
  
  <pre><code> (define (square x) (* x x))
  
   (define remainder mod)
  
   (define (miller-rabin-expmod base exp m) 
     (define (squaremod-with-check x) 
       (define (check-nontrivial-sqrt1 x square) 
         (if (and (= square 1) 
                  (not (= x 1)) 
                  (not (= x (- m 1)))) 
             0 
             square)) 
       (check-nontrivial-sqrt1 x (remainder (square x) m))) 
     (cond ((= exp 0) 1) 
           ((even? exp) (squaremod-with-check 
                         (miller-rabin-expmod base (/ exp 2) m))) 
           (else 
            (remainder (* base (miller-rabin-expmod base (- exp 1) m)) 
                       m))))
  
  (define expmod miller-rabin-expmod)
  
  (print (expmod 1 14 15))
  (print (expmod 2 14 15))
  (print (expmod 3 14 15))
  (print (expmod 4 14 15))
  (print (expmod 5 14 15))
  (print (expmod 6 14 15))
  (print (expmod 7 14 15))
  (print (expmod 8 14 15))
  (print (expmod 9 14 15))
  (print (expmod 10 14 15))
  (print (expmod 11 14 15))
  (print (expmod 12 14 15))
  (print (expmod 13 14 15))
  (print (expmod 14 14 15))
  </code></pre>
  
  <p>Edit 2: I have now also manually calculated the steps of <code>expmod(a, 14, 15)</code> (which always recurses via <code>exp = 14</code>, <code>exp = 7</code>, <code>exp = 6</code>, <code>exp = 3</code>, <code>exp = 2</code>, <code>exp = 1</code>, <code>exp = 0</code>), for all values of <code>a</code> from 1 to 14, and I'm certain that only <code>a = 4</code> and <code>a = 11</code> encounter a nontrivial square root of 1. So I'm inclined to think that SICP is either wrong about this, or is not expressing itself clearly.</p>