## 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 *a ^{d}* ≡ 1 (mod

*n*) and the inner

`loop`

indexes through the factors of 2 testing the other condition *a*

^{d·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.

This is where the criticism of Python being too slow gains ground, in my

opinion. Trying to do number intensive stuff takes a while. My solution

has a few changes: I went with probabilistic testing instead of deterministic

testing, had my

`_factors(n)`

give just primes, etc. The logic ismostly the same.