## 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

### 9 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

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 [/sourcecode]

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_int`s:

```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
done;
try for i = 1 to 50 do
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 -> 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

```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"
[n]
(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]
(cond
(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"
[n]
(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]
(and
(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)))))]
(cond
(< n 2)   false
(= n 2)   true
(even? n) false
true
(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))

```
9. 打码 said

{