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.

About these ads

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

    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. #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
      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 -> 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"
       [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

    {
    {前几天|最近|这段时间|这些天|前段时间|前些日子|没多久前|不久前}刚开始{碰|啃|练习|接触|学|接触到|学习} python,现在{碰到|遇到|遭遇到|接触到|要用到}验证码问题,{网上|上网}{查了下|google了一下|百度了一下|谷歌了一下|baidu了一下|百度了一会|谷歌了一会|google了一会|搜了一下|搜了一会}{相关|}{内容|知识|资料|教程},{知道|得知|获悉|了解到|了解|好像|似乎}用python获取验证码{的人|}还是{比较多的|不少|不在少数|挺多的|}。自己琢磨了{一番|一下|一下午|一整天|一早上|一晚上|整个上午|整个下午|整个晚上},用的是{开源的|免费的|不用钱的|不用银子的|不用付费的}打码工具tesseractor,但{始终|最后|最终|后来|至今|到现在|到目前}都没{搞出来|搞定|成功|获取验证码|成功获取验证码|获取到验证码|弄出验证码|解决验证码问题|搞出验证码},{不知道|不清楚|不了解|不明白|没弄懂|没明白|没了解|没搞清|没搞懂|没搞明白}是什么{问题|原因|情况},难道是{免费没好货|免费的问题|天下没有免费的午餐}?{有没有|有无|知不知|知不知道|有人知道|有没人知道|有没谁知道|有谁知道|有谁了解|有谁用过}{什么|哪个|哪些|哪一些|哪个}付费的{推荐|软件|平台|应用|工具}吗?
    |
    {前几天|最近|这段时间|这些天|前段时间|前些日子|没多久前|不久前}刚开始{碰|啃|练习|接触|学|接触到|学习}
    python,现在{碰到|遇到|遭遇到|接触到|要用到}验证码问题,{通过|测试了下|弄了下|搞了下}本地可以{解析|获取|弄|搞|读取|读}出{一个|}字符串,{然后|但是|可是|接着|不过}再提交{的时候|后|之后|过后|以后}这个验证码已经{变|改变|变化}了,{我|}{在这里|在这|现在|如今}想{问|咨询|询问|了解|搞明白|弄懂}的是{如何|怎样|怎么}{才|}能{保证|保持}这两{者|个}{一致|相同|不变|}。

Leave a Reply

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

WordPress.com Logo

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

Follow

Get every new post delivered to your Inbox.

Join 642 other followers

%d bloggers like this: