Proving Primality

February 2, 2010

[ Updated 14 April 2010 to fix a bug. ]

We use the trial-division function td-factors from the exercise on wheel factorization and expm from the Standard Prelude:

(define (prime? n)
  (if (even? n) (= n 2)
    (let* ((n-1 (- n 1)) (fs (td-factors n-1)))
      (let loop1 ((b 2))
        (cond ((= b n) #f)
              ((= (expm b n-1 n) 1)
                (let loop2 ((qs fs))
                  (cond ((null? qs) #t)
                        ((= (expm b (/ n-1 (car qs)) n) 1)
                          (loop1 (+ b 1)))
                        (else (loop2 (cdr qs))))))
              (else (loop1 (+ b 1))))))))

As an example, we compute the prime factors of 289-1. The prime factors of n-1 are 2, 3, 5, 17, 23, 89, 353, 397, 683, 2113, and 2931542417, and it is 2(n-1)/89 ≡ 4096 (mod n) that provides the counterexample that proves 289-1 is prime:

> (prime? (- (expt 2 89) 1))
#t

Since the first version of this program had a bug, we provide the following test program, which uses the primes function from a previous exercise:

(define (test-prime? n)
  (let ((ps (primes n)))
    (do ((i 2 (+ i 1))) ((= i n))
      (let ((p? (prime? i)))
        (when (and p? (not (member i ps)))
          (display i) (display " found to be prime") (newline))
        (when (and (not p?) (member i ps))
          (display i) (display " found to be composite") (newline))))))

Calling (test-prime? 1000) shows no errors.

The speed of this algorithm depends primarily on the speed by which n-1 can be factored, so it is better to replace the factorization by trial division with one of the other factorization functions that we have studied. But if you do that, be sure that if you use one of the algorithms that can return a non-prime factor, you don’t use a probabilistic prime checker to determine whether or not the factor is prime!

You can run this program at http://programmingpraxis.codepad.org/1ttTLsZF.

About these ads

Pages: 1 2

6 Responses to “Proving Primality”

  1. [...] Praxis – Proving Primality By Remco Niemeijer In today’s Programming Praxis exercise we have to implement an algorithm to prove the primality of a number. [...]

  2. Remco Niemeijer said

    My Haskell solution (see http://bonsaicode.wordpress.com/2010/02/02/programming-praxis-proving-primality/ for a version with comments):

    import Data.Bits
    import Data.List
    
    tdFactors :: Integer -> [Integer]
    tdFactors n = f n [2..floor . sqrt $ fromIntegral n] where
        f _ []     = []
        f r (x:xs) | mod r x == 0 = x : f (div r x) (x:xs)
                   | otherwise    = f r xs
    
    expm :: Integer -> Integer -> Integer -> Integer
    expm b e m = foldl' (\r (b', _) -> mod (r * b') m) 1 .
                 filter (flip testBit 0 . snd) .
                 zip (iterate (flip mod m . (^ 2)) b) .
                 takeWhile (> 0) $ iterate (`shiftR` 1) e
    
    isPrime :: Integer -> Bool
    isPrime n = f 2 $ tdFactors (n - 1) where
        f _ []     = False
        f b (q:qs) | expm b (n - 1)         n /= 1 = f (b + 1) (q:qs)
                   | expm b (n - 1 `div` q) n /= 1 = True
                   | otherwise                     = f b qs
    
  3. Mike said

    I found the description of the algorithm confusing. When I implemented it, the routine would report some multiples of 5 (e.g., 15 and 35) to be prime.

    Studying the proof at the mathpages.com link, given above, I think the algorithm should be to find a b so that b^(x-1) == 1 (mod n) then, if b^((n-1)/q) != 1 (mod n) for *all* prime divisors q of n-1, n is prime. If b^((n-1)/q) == 1 (mod n) for some q, then try a new value for b. This matches the pseudocode at wikipedia.

    My python code is below:

    from itertools import chain, cycle
    
    def factors_of(n):
        f = 2
        for step in chain([0,1,2,2], cycle([4,2,4,2,4,6,2,6])):
            f += step
            
            if f*f > n:
                if n != 1:
                    yield n
                break
            
            if n%f == 0:
                yield f
                while n%f == 0:
                    n /= f
    
                    
    def is_prime(n):
        '''proves primality of n using Lucas test.'''
        x = n-1
    
        if n == 2 or n%2 == 0:
            return False
    
        factors = list( factors_of(x))
        
        b = 2
        while b < n:
            if pow( b, x, n ) == 1:
    
                for q in factors:
                    if pow( b, x/q, n ) == 1:
                        break
                else:
                    return True
                
            b += 1
            
        return False
    
  4. programmingpraxis said

    It wasn’t confusing, it was wrong. Thanks for pointing out the bug. I’ve updated the solution and added a test function.

    By the way, your version of is_prime incorrectly reports that 2 is composite. You should change the test at lines 22 and 23 to read:

    if n%2 == 0:
        return n == 2

    Don’t feel bad; Jon Bentley got that wrong in one of his books, so you’re in good company.

  5. Mike said

    D’oh…guess that’ll teach me not to edit a routine and not test it again before posting.
    Is there a way to edit my earlier post fix lines 22 and 23? Or should I just repost a correct version?

    My confusion was over the word “any” in

    “… if b^((n-1)/q) != 1 (mod n) for *any* prime divisor q of n-1, n is prime.”

    At first, I read the phrase to mean if I found any one q for which the non-congruence was true, then n is prime. However, the correct reading is that the non-congruence holds no matter which q you try. Therefore, you have to try every q to prove n is prime. Perhaps using “every” instead of “any” would be clearer.

  6. Mike said

    Here’s a corrected version of the python code for proving a number is prime using the Lucas algorithm.

    
    from itertools import chain, cycle
    
    def factors_of(n):
        f = 2
        for step in chain([0,1,2,2], cycle([4,2,4,2,4,6,2,6])):
            f += step
    
            if f*f > n:
                if n != 1:
                    yield n
                break
    
            if n%f == 0:
                yield f
                while n%f == 0:
                    n /= f
    
    def is_prime(n):
        '''proves primality of n using Lucas test.'''
        x = n-1
    
        if n%2 == 0:
            return n == 2
    
        factors = list( factors_of(x))
    
        b = 2
        while b < n:
            if pow( b, x, n ) == 1:
    
                for q in factors:
                    if pow( b, x/q, n ) == 1:
                        break
                else:
                    return True
    
            b += 1
    
        return False
    

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 621 other followers

%d bloggers like this: