SICP Exercise 1.27: Carmichael Numbers

Demonstrate that the Carmichael numbers listed in footnote 47 really do fool the Fermat test. That is, write a procedure that takes an integer n and tests whether aⁿ is congruent to a modulo n for every a<n, and try your procedure on the given Carmichael numbers.

Footnote 47 states that:

"Numbers that fool the Fermat test are called Carmichael numbers, and little is known about them other than that they are extremely rare. There are 255 Carmichael numbers below 100,000,000. The smallest few are 561, 1105, 1729, 2465, 2821, and 6601. In testing primality of very large numbers chosen at random, the chance of stumbling upon a value that fools the Fermat test is less than the chance that cosmic radiation will cause the computer to make an error in carrying out a "correct" algorithm. Considering an algorithm to be inadequate for the first reason but not for the second illustrates the difference between mathematics and engineering.

So we want to check that the following numbers fool the Fermat test: 561, 1105, 1729, 2465, 2821, and 6601.

Before we do that, we need a procedure that performs the Fermat test for all values in the range 1≤a<n. First, let's look at the original fermat-test implementation:

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

This has an internal procedure definition for try-it, which tests whether or not the Fermat test holds for a single integer, a, which fermat-test calls using a random value in the range 1≤a<n. We can use this as the basis for a new procedure (which I've called it fermat-test-all) by removing the randomness and replacing it with iterative testing of all of the values in the range 1≤a<n. To do this we:

• Modify try-it that:
  1. It checks to see if the Fermat test holds for the current value of a.
  2. If it does then it calls itself for the next value of a.
  3. If it doesn't then it returns #f.
  4. If it reaches a=n then it returns #t, as the Fermat test has succeeded for all tested numbers.
• Modify the outer call to try-it so that it starts at 1.

Here's the updated procedure:

(define (fermat-test-all n)
  (define (try-it a)
    (cond ((= a n) #t)
          ((= (expmod a n n) a) (try-it (+ a 1)))
          (else #f)))
  (try-it 1))

I also produced a procedure charmichael? that checks whether or not a given number is a Charmichael number or not. I.e. it checks to see that it passes the Fermat test and it checks that the number is not prime:

(define (charmichael? n)
  (and (fermat-test-all n) (not (prime? n))))

When we run this for the given Charmichael numbers we get:

> (charmichael? 561)
#t
> (charmichael? 1105)
#t
> (charmichael? 1729)
#t
> (charmichael? 2465)
#t
> (charmichael? 2821)
#t
> (charmichael? 6601)
#t

...and just as a quick sanity check:

> (charmichael? 1)
#f
> (charmichael? 2)
#f
> (charmichael? 3)
#f
> (charmichael? 4)
#f
> (charmichael? 5)
#f
> (charmichael? 6)
#f
> (charmichael? 7)
#f

SICP Exercise 1.26: A Broken expmod

Louis Reasoner is having great difficulty doing exercise 1.24. His fast-prime? test seems to run more slowly than his prime? test. Louis calls his friend Eva Lu Ator over to help. When they examine Louis's code, they find that he has rewritten the expmod procedure to use an explicit multiplication, rather than calling square:

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (* (expmod base (/ exp 2) m)
                       (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

"I don't see what difference that could make," says Louis. "I do." says Eva. "By writing the procedure like that, you have transformed the Θ(log n) process into a Θ(n) process." Explain.

The technique of using successive squaring (introduced in section 1.2.4) to reduce an Θ(n) process to a Θ(log n) process is designed to reduce the number of calculations performed. This technique is applied in expmod so that it calculates (expmod base (/ exp 2) m) only once and then squares the result. The equivalent (Θ(n)) approach would be to calculate (expmod base (- exp 1) m) in each recursive call, regardless of whether exp was even or not.

The change made by Louis causes the process to calculate (expmod base (/ exp 2) m) twice instead of once for each invocation of (expmod base exp m), and then multiplies the results together. Note that I said "for each invocation of (expmod base exp m)" quite deliberately here...

In Louis' version, in order to calculate (expmod base exp m) the interpreter has to first calculate (expmod base (/ exp 2) m) twice. However, for each calculation of (expmod base (/ exp 2) m) the interpreter first has to calculate (expmod base (/ exp 4) m) twice, making four calculations of (expmod base (/ exp 4) m), and so on. The recursion is still log n deep, but the number of calculations at level i in the recursion is 2ⁱ. Here's how the number of calls to expmod pan out for the different processes:

exp	Θ(log n) Process	Θ(n) Process	Louis' Process
0	1	1	1
1	2	2	3
2	3	3	7
4	4	5	15
8	5	9	31
16	6	17	63
32	7	33	127

So we can see that the Θ(log n) process takes log₂n steps, and the Θ(n) process takes n + 1 steps. What about Louis' process? Well, that takes 4n - 1 steps, which is Θ(n). Note that the original prime? test is Θ(√n), so Louis is right - his fast-prime? test does run more slowly than his prime? test.

SICP Exercise 1.25: A Simpler expmod?

Alyssa P. Hacker complains that we went to a lot of extra work in writing expmod. After all, she says, since we already know how to compute exponentials, we could have simply written

(define (expmod base exp m)
  (remainder (fast-expt base exp) m))

Is she correct? Would this procedure serve as well for our fast prime tester? Explain.

Yes. And no.

Yes, we could have written what Alyssa P. Hacker suggests and it would have successfully performed the prime test. But no, it wouldn't serve as well for our fast prime tester - it would have taken much longer to perform the test.

Why? Well, as I discussed in my solution to exercise 1.24, once the numbers exceed a certain value (2³⁰ on my interpreter) the interpreter needs to use a different representation that can handle larger values. And when it does this, the performance of the calculations degrades. In this case we're raising numbers to powers over 1,000, 10,000, etc, and so we'll quickly reach the point where this change occurs.

Here's what happened when I plugged it into my code from exercise 1.24:

> (timed-prime-test 1009)
1009 *** 47611.0
> (timed-prime-test 1013)
1013 *** 49780.0
> (timed-prime-test 1019)
1019 *** 49577.0
> (timed-prime-test 10007)
10007

...and I gave up waiting after several minutes. We've gone from approximately 1.6 seconds to around 48 seconds each of the first three primes over 1,000!

It's possibly worth noting that one of the footnotes in section 1.2.6 pretty much points out the flaw in using fast-expt:

"The reduction steps in the cases where the exponent e is greater than 1 are based on the fact that, for any integers x, y, and m, we can find the remainder of x times y modulo m by computing separately the remainders of x modulo m and y modulo m, multiplying these, and then taking the remainder of the result modulo m. For instance, in the case where e is even, we compute the remainder of b^e/2 modulo m, square this, and take the remainder modulo m. This technique is useful because it means we can perform our computation without ever having to deal with numbers much larger than m."

The magic phrase being "it means we can perform our computation without ever having to deal with numbers much larger than m". When we don't use this technique we do end up dealing with numbers larger than m. Much larger. The interpreter then ends up using less performant integer representations for most of its calculations and we get the degredation witnessed above.