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
[...] 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))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.
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_intbetween 0 and a specifiedmax(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 eFinally, 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 done; 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 done; if not !any then raise Exit done; true with Exit -> falseTwo 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.
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);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