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.

About these ads

Pages: 1 2

One Response to “Carmichael Numbers”

  1. Graham said

    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 is
    mostly the same.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 576 other followers

%d bloggers like this: