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 ≤ an-1 and check if as ≡ 1 (mod n) or a2j s ≡ -1 (mod n) for some 0 ≤ jr-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.

About these ads

Pages: 1 2

8 Responses to “Primality Checking”

  1. [...] 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

    My Haskell solution (see 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. #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))
  4. Graham said

    I’m slowly working my way through old exercises now that I have some free time. Here’s
    my attempt in Common Lisp; I’m trying to learn the language, even though I already had Pythoned Miller-Rabin previously.

  5. Wow, so much code needed in OCaml to solve this exercise in a self-contained fashion! First of all, a general-purpose function to return a random Big_int between 0 and a specified max (exclusive):

    let random_big_int =
      let open Big_int in
      let open Nat in
      let random_limb = match length_of_digit with
      | 64 -> fun nat ofs tmp ->
        set_digit_nat nat ofs (Random.bits ());
        shift_left_nat nat ofs 1 tmp 0 30;
        set_digit_nat tmp 0 (Random.bits ());
        lor_digit_nat nat ofs tmp 0;
        shift_left_nat nat ofs 1 tmp 0 30;
        set_digit_nat tmp 0 (Random.bits () land 7);
        lor_digit_nat nat ofs tmp 0
      | 32 -> fun nat ofs tmp ->
        set_digit_nat nat ofs (Random.bits ());
        shift_left_nat nat ofs 1 tmp 0 30;
        set_digit_nat tmp 0 (Random.bits () land 3);
        lor_digit_nat nat ofs tmp 0
      | _  -> assert false
      in fun max ->
      let nat = nat_of_big_int max in
      let len = num_digits_nat nat 0 (length_nat nat) in
      let res = create_nat len
      and tmp = create_nat 1 in
      for i = 0 to len - 1 do random_limb res i tmp done;
      mod_big_int (big_int_of_nat res) max

    (You can’t go much lower level than this.) Next, some utility functions on Big_ints:

    let is_zero_big_int n = Big_int.sign_big_int n == 0
    let is_even_big_int n = Big_int.( is_zero_big_int (and_big_int n unit_big_int) )
    let modsquare_big_int x n = Big_int.( mod_big_int (square_big_int x) n)
    let modexp_big_int x e n =
      let open Big_int in
      let rec go y z e =
        if is_zero_big_int e then y else
        let y = if is_even_big_int e
          then y
          else mod_big_int (mult_big_int y z) n
        in go y (modsquare_big_int z n) (shift_right_big_int e 1)
      in go unit_big_int x e

    Finally, an imperative-style Rabin-Miller test:

    let is_prime n =
      let open Big_int in
      if le_big_int n unit_big_int || is_even_big_int n
        then invalid_arg "is_prime" else
      let m = pred_big_int n in
      let r = ref 0
      and s = ref m in
      while is_even_big_int !s do
        incr r;
        s := shift_right_big_int !s 1
      try for i = 1 to 50 do
        let a = add_int_big_int 2 (random_big_int (add_int_big_int (-2) n)) in
        let x = ref (modexp_big_int a !s n)
        and j = ref !r in
        let any = ref (eq_big_int !x unit_big_int) in
        while !j != 0 && not !any do
          if eq_big_int !x m
            then any := true
            else x := modsquare_big_int !x n;
          decr j
        if not !any then raise Exit
      with Exit -> false

    Two optimizations were applied: first, the random a should be greater than 1; second, the successive powers of a in the inner loop are calculated inductively by (modular) squaring.

  6. I made a mistake in random_big_int. Lines 10 and 11 are wrong, they should read:

        shift_left_nat nat ofs 1 tmp 0 4;
        set_digit_nat tmp 0 (Random.bits () land 15);
  7. David said

    This is a Haskell program I wrote (ironically enough around May 2009, though I hadn’t heard about this website back then.) This was my first relatively non-trivial Haskell program.

    import System.Random
    import System (getArgs)
    -- return number N as (s,d) where N = (2^S)*d
    pow2factor n
       | odd n  = (0, n)
       | even n = (s+1, d) where (s,d) = pow2factor (n `div` 2)
    -- compute (a^n) (mod m)
    --     (efficient algorithm, doing modulo arithmetic after each step)
    square x = x*x
    pow a 0 m = 1
    pow a n m
       | odd  n = (a * (pow a (n - 1) m)) `mod` m
       | even n = (square (pow a (n `div` 2) m)) `mod` m
    -- Miller-Rabin: given an integer n, determine if it is prime with
    -- probability 4 ^ -k, where k is the number of tests performed.
    probablyPrime n tests =
        let (s, d) = pow2factor (n-1) in not (any (\x -> testComposite x s d) tests) where
    	testComposite a s d =
    	    (pow a d n) `notElem` [1, n-1] &&
    	    all (\r -> (pow a (r*d) n) /= n-1) (take (s-1) (iterate (2*) 2))
    -- main
    testPrime n =
        do rnd <- newStdGen
           return (probablyPrime n (take 10 (randomRs (2, n-2) rnd)))
    genprime :: Integer -> IO Integer
    genprime bits =
        do rnd <- getStdGen
           thePrime <- firstprime (randomRs (2, 2^(bits-1)) rnd)
           return thePrime where
               firstprime (h:t) =
    	       let candidate = 2*h+1 in
    		   do isPrime <- testPrime candidate
    		      if isPrime then return candidate else firstprime t
    main =
       do args <- getArgs
          case args of
    	 ["-g", bits] -> do x <- genprime ((read bits) :: Integer)
    	                    putStrLn (show x)
    	 ["-t", n] -> do isPrime <- testPrime ((read n) :: Integer)
    			 putStrLn (if isPrime then "prime" else "composite")
    	 otherwise -> error "Usage: miller-rabin [-g bits] [-t number]"

    To test 2^89 -1 use,
    .\miller-rabin.exe -t 618970019642690137449562111

  8. David said

    Rewrote the Haskell Miller-Rabin in Clojure as my first Clojure program. It is line-by-line translation, pretty much. Main difference is that one can’t generate random numbers > 2^63, so that is taken into account when generating tests.

    (defn factor-2s
       "return argument N as vector [s,d] where N = (2^S)*d"
       (loop [n' n  e 0]
          (if (even? n') 
             (recur (/ n' 2) (inc e))
             [e n'])))
    (defn square [x]  (* x x))
    (defn pow
       "compute (a^n) (mod m)
           (efficient algorithm, doing modulo arithmetic after each step)"
       [a n m]
          (zero? n)  1
          (odd? n)   (mod (* a (pow a (dec n) m)) m)
          (even? n)  (mod (square (pow a (/ n 2) m)) m)))
    (defn random-gen
       "Lazily computed list of random numbers"
       (map #(long (rand %1)) (cycle [n])))
    (defn probably-prime
       "Miller-Rabin: given an integer n, determine if it is prime with
        probability 1 - (4^-k), where k is the number of tests performed."
       [accuracy n]
       (let [[s d] (factor-2s (dec n))
            composite-witness?  (fn [a]
                  (let [t (pow a d n)]  (and (not= t 1) (not= t (dec n))))
                  (every? (fn [r] (not= (pow a (* r d) n) (dec n)))
                          (take (dec s) (iterate #(+ %1 %1) 2)))))]
             (< n 2)   false
             (= n 2)   true
             (even? n) false
                (not-any? composite-witness?
                      (take accuracy (map (partial + 2) (random-gen (min (- n 4) 0x7FFFFFFFFFFFFFFF))))))))
    (def prime-accuracy 10)
    (def prime?  (partial probably-prime prime-accuracy))

Leave a Reply

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

You are commenting using your 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


Get every new post delivered to your Inbox.

Join 576 other followers

%d bloggers like this: