Proving Primality

February 2, 2010

We examined algorithms for checking the primality of a number in two previous exercises. Both algorithms are probabilistic; the Rabin-Miller algorithm has known false positives, and Carl Pomerance proved that the Baillie-Wagstaff algorithm has an infinite number of false positives, even though none are known.

Today’s exercise describes an algorithm for proving that a number is prime, not probably but certainly. The algorithm works by a theorem of Edouard Lucas, which is based on Fermat’s Little Theorem (if n is prime then bn-1 ≡ 1 (mod n) for every integer b coprime to n):

If, for some integer b, bn-1 ≡ 1 (mod n), and if b(n-1)/q ≢ 1 (mod n) for any prime divisor q of n-1, then n is a prime number.

Thus, if the factors of n-1 are known, it is easy to determine the primality of n.

Your task is to write a function that determines the primality of an integer using Lucas’ primality test; use it to prove the primality of 289-1. When you are finished, you are welcome to read or run a suggested solution , or to post your own solution or discuss the exercise in the comments below.

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

%d bloggers like this: