Carmichael Numbers

December 28, 2010

We begin with Fermat’s pseudo-prime test to a single base, which follows the formula exactly:

```(define (fermat-pseudo-prime? a n)   (and (= (gcd a n) 1)        (= (expm a (- n 1) n) 1)))```

Here’s an example. Note that 561 is not a Fermat pseudo-prime to 3 because 3 is a factor of 561:

```> (fermat-pseudo-prime? 2 561) #t > (fermat-pseudo-prime? 3 561) #f```

The strong pseudo-prime test is somewhat more complicated. The outer `loop` extracts factors of 2 from n−1, the second `if` tests the condition that ad ≡ 1 (mod n) and the inner `loop` indexes through the factors of 2 testing the other condition ad·2r ≡ −1 (mod n):

```(define (strong-pseudo-prime? a n)   (let loop ((r 0) (s (- n 1)))     (if (even? s) (loop (+ r 1) (/ s 2))       (if (= (expm a s n) 1) #t         (let loop ((j 0) (s s))           (cond ((= j r) #f)                 ((= (expm a s n) (- n 1)) #t)                 (else (loop (+ j 1) (* s 2)))))))))```

Here are the same two examples above; note that 2 is now a witness to the compositeness of 561, whereas it failed with the Fermat pseudo-prime test:

```> (strong-pseudo-prime? 2 561) #f > (strong-pseudo-prime? 3 561) #f```

Either test can be turned into a primality test, though of course the strong pseudo-prime test is preferable to the Fermat pseudo-prime test. We could use random bases, as in the Miller-Rabin test, but instead we choose a deterministic test using the primes less than 200:

```(define pseudo-prime?   (let ((as (primes 200)))     (lambda (method? n)       (if (member n as) #t         (let loop ((as as))           (cond ((null? as) #t)                 ((not (method? (car as) n)) #f)                 (else (loop (cdr as)))))))))```

Given that higher-order function, the two primality tests are easy:

```(define (fermat-prime? n)   (pseudo-prime? fermat-pseudo-prime? n))```

```(define (strong-prime? n)   (pseudo-prime? strong-pseudo-prime? n))```

Here are some examples:

```> (fermat-prime? 561) #f > (strong-prime? 561) #f```

We come now to the two Carmichael tests. First, using Fermat pseudo-primes:

```(define (carmichael? n)   (let ((ps (primes n)))     (if (= (car (reverse ps)) n) #f       (let loop ((ps (primes n)))         (cond ((null? ps) #t)               ((not (= (gcd (car ps) n) 1)) (loop (cdr ps)))               ((not (fermat-pseudo-prime? (car ps) n)) #f)               (else (loop (cdr ps))))))))```

Note that this differs from the Fermat pseudo-prime test by first checking that n is coprime to the base, which is part of the definition of Carmichael numbers. We also check that n is prime by sieving, because a number is only a Carmichael number if it is composite, and of course the Fermat pseudo-prime test can’t make that distinction. Here’s the list of Carmichael numbers less than 10000:

```> (let loop ((c 10000) (cs '()))     (cond ((= c 2) cs)           ((carmichael? c)             (loop (- c 1) (cons c cs)))           (else (loop (- c 1) cs)))) (561 1105 1729 2465 2821 6601 8911)```

The Carmichael test based on Korselt’s criterion relies on the factorization of n; the `(pair? (cdr fs))` test discards primes that have only themselves as factors, and the `(apply < fs)` fails if any factors are repeated, meaning that n is not square-free.

```(define (carmichael? n)   (let ((fs (factors n)))     (and (odd? n) (pair? (cdr fs)) (apply < fs)          (let loop ((fs fs))            (cond ((null? fs) #t)                  ((not (zero? (modulo (- n 1) (- (car fs) 1)))) #f)                  (else (loop (cdr fs))))))))```

We need a function to factor integers, so we choose this simple implementation of trial division, knowing that it can be replaced with a more sophisticated method if needed:

```(define (factors n)   (if (even? n) (cons 2 (factors (/ n 2)))     (let loop ((n n) (f 3) (fs '()))       (cond ((< n (* f f)) (reverse (cons n fs)))             ((zero? (modulo n f)) (loop (/ n f) f (cons f fs)))             (else (loop n (+ f 2) fs))))))```

The list of Carmichael numbers is the same as before:

```> (let loop ((c 10000) (cs '()))     (cond ((= c 2) cs)           ((carmichael? c)             (loop (- c 1) (cons c cs)))           (else (loop (- c 1) cs)))) (561 1105 1729 2465 2821 6601 8911)```

We used `expm` from the Standard Prelude and `primes` from the exercise on the Sieve of Eratosthenes. You can run the code at http://programmingpraxis.codepad.org/UudIdev3. Sloane’s A002997 gives the Carmichael numbers.

testing, had my `_factors(n)` give just primes, etc. The logic is