Primality Checking
May 1, 2009
In a previous exercise you wrote a function that returned a list of prime numbers, and in another exercise you used that function to find a particular prime number. This exercise looks at prime numbers from a different perspective by considering a function that takes a number and determines if it is prime.
The algorithm that we will consider was developed initially by Gary Miller and refined by Michael Rabin, and is probabilistic in nature. It works like this: Express the odd number n to be factored as n = 2r s + 1 with s odd. Then choose a random integer a with 1 ≤ a ≤ n-1 and check if as ≡ 1 (mod n) or a2j s ≡ -1 (mod n) for some 0 ≤ j ≤ r-1. (Some browsers render that last equation poorly; it’s a raised to the power 2 to the j times s.) A prime number will pass the check for all a. A composite number will pass the check for about 1/4 of the possible as and fail the check for the remaining 3/4 of the possible as. Thus, to determine if a number n is prime, check multiple as; if k as are checked, this algorithm will err less than one time in 4k. Most primality checkers set k to somewhere between 25 and 50, making the chance of error very small.
Your task is to write a function that determines if an input number n is prime, then to determine if 289-1 is prime. When you are finished, you are welcome to read or run a suggested solution, or post your solution or discuss the exercise in the comments below.
Pages: 1 2
[...] Praxis – Primality Checking By Remco Niemeijer Today’s Programming Praxis problem is about checking whether or not a number is prime. We’re supposed [...]
My Haskell solution (see http://bonsaicode.wordpress.com/2009/05/01/programming-praxis-primality-checking/ for a version with comments):
import Control.Arrow import Data.Bits import Data.List import System.Random isPrime :: Integer -> StdGen -> Bool isPrime n g = let (s, d) = (length *** head) . span even $ iterate (flip div 2) (n - 1) xs = map (expm n d) . take 50 $ randomRs (2, n - 2) g in all (\x -> elem x [1, n - 1] || any (== n - 1) (take s $ iterate (expm n 2) x)) xs expm :: Integer -> Integer -> Integer -> Integer expm m e b = foldl' (\r (b', _) -> mod (r * b') m) 1 . filter (flip testBit 0 . snd) . zip (iterate (flip mod m . (^ 2)) b) $ takeWhile (> 0) $ iterate (flip shiftR 1) e main :: IO () main = print . isPrime (2 ^ 89 - 1) =<< getStdGen#lang scheme (require srfi/1 ; list-tabulate srfi/27) ; random-integer (define (factor2 n) ; return r and s s.t n = 2^r * s where s is odd (let loop ([r 0] [s n]) ; invariant: n = 2^r * s (let-values ([(q r) (quotient/remainder s 2)]) (if (zero? r) (loop (+ r 1) q) (values r s))))) (define (miller-rabin n) ; Input: n odd Output: n prime? (define (mod x) (modulo x n)) (define (^ x m) (cond [(zero? m) 1] [(even? m) (mod (sqr (^ x (/ m 2))))] [(odd? m) (mod (* x (^ x (- m 1))))])) (define (check? a) (let-values ([(r s) (factor2 (sub1 n))]) (and (member (^ a s) (list 1 (mod -1))) #t))) (andmap check? (list-tabulate 50 (λ (_) (+ 2 (random-integer (- n 3))))))) (define (prime? n) (cond [(< n 2) #f] [(= n 2) #t] [(even? n) #f] [else (miller-rabin n)])) (prime? (- (expt 2 89) 1))