Primality Checking

May 1, 2009

Although the math looks difficult, the code is quite straight forward. Check? calculates r and s in the outer loop, checks if as is 1 modulo n, and then the inner loop checks the other condition for each j from 0 to r-1.

(define (check? 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)))))))))

Then prime? tests a few boundary conditions and performs fifty random calls to check?.

(define (prime? n)
  (cond ((< n 2) #f) ((= n 2) #t) ((even? n) #f)
        (else (let loop ((k 50))
                (cond ((zero? k) #t)
                      ((not (check? (randint 1 n) n)) #f)
                      (else (loop (- k 1))))))))

Expm (modular exponentiation) and randint come from the Standard Prelude. A quick test shows that 289-1 = 618970019642690137449562111 is prime:

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

You can run this code at http://codepad.org/rSDxFrZn.

Pages: 1 2

Posted by programmingpraxis
Filed in Exercises
3 Comments »

3 Responses to “Primality Checking”

  1. Programming Praxis - Primality Checking « Bonsai Code said
    May 1, 2009 at 2:19 PM

    [...] Praxis – Primality Checking By Remco Niemeijer Today’s Programming Praxis problem is about checking whether or not a number is prime. We’re supposed [...]

  2. Remco Niemeijer said
    May 1, 2009 at 2:20 PM

    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
  3. Jens Axel Søgaard said
    May 4, 2009 at 6:49 PM 
    #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))

Leave a Reply