# 快速素数判定

Article Outline

``````(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))

;return the smallest divisor of n
(define (smallest-divisor n)
;find a divisor(for 2 to n) of n
(define (find-divisor n test-divisor)
(cond ((< n (sqrt test-divisor)) n)
((= (remainder n test-divisor) 0) test-divisor)
(else (find-divisor n (+ test-divisor 1)))))

(find-divisor n 2))

;check n is prime or not
(define (prime? n)
(= n (smallest-divisor n )))``````

``````a^n = (a^(n/2))^2; 当n为偶数
a^n = a* a^(n-1); 当n为奇数``````

``````a^n mod m = ((a^(n / 2))^2) mod m;  当n为偶数
a^n mod m = (a* a^(n-1)) mod m;     当n为奇数``````

``````(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))

(define (expmod base exps m)
(define (square n) (* n n))
(cond ((= exps 0) 1)
((even? exps)
(remainder (square (expmod base (/ exps 2) m))
m))
(else
(remainder (* base (expmod base (- exps 1) m))
m))))

;fermat test, check n. if (a^n) mod n == a then in all probability n is a prime
(define (fermat-test n)
(define (try-it a)
(= a (expmod a n n)))
;(= a (remainder (expt a n) n)))
(try-it (+ 1 (random (- n 1)))))

;get prime by fermat-test. test n is prime or not, try times time
(define (fermat-prime? n times)
(cond ((= 0 times) #t)
((fermat-test n) (fermat-prime? n (- times 1)))
(else #f)))``````

561，1105，1729，2465，2821，6601，8911

Matrix67：Miller-Rabin 素性测试同样是不确定算法，我们把可以通过以a为底的 Miller-Rabin 测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以 2 为底的强伪素数为 2047。第一个以 2 和 3 为底的强伪素数则大到 1 373 653。

``````(define (remainder num divisor)
(- num (* divisor (round (/ num divisor)))))

(define (mr-expmod base exps m)
(define (square x) (* x x))
(define (check-square-root x n)
(define (check r)
(if (and (not (= x 1)) (not (= x (- n 1))) (= r 1)) 0 r))
(check (remainder (square x) n)))

(cond ((= exps 0) 1)
((even? exps)
(check-square-root(mr-expmod base (/ exps 2) m) m))
(else
(remainder (* base (mr-expmod base (- exps 1) m))
m))))

;Miller-Rabin test, check n. if (a^(n-1)) mod n == 1 then in all probability n is a prime
(define (mr-test n)
(define (try-it a)
(= 1 (mr-expmod a (- n 1) n)))
(try-it (+ 1 (random (- n 1)))))

;get prime by Miller-Rabin-test. test n is prime or not, try times time
(define (mr-prime? n times)
(cond ((= 0 times) #t)
((mr-test n) (mr-prime? n (- times 1)))
(else #f)))
``````

``````    常规方法           常规方法的因数2优化版本       Miller-Rabin测试法
[1000, 1100]           1 ms                1 ms                     1 ms
[10000, 10100]         13 ms               8 ms                     1 ms
[100000, 100100]       142 ms              78 ms                    1 ms
[1000000, 1000100]     1337 ms             830 ms                   1 ms``````

``````    \$ ./prime
Normal:  return 0 cost 10ms
Fermat: return 0 cost 24ms
Miller-Rabin: return 0 cost 28ms``````

``````    \$ ./prime
Normal: return 0 cost 3088ms
Fermat: return 0 cost 824ms
Miller-Rabin: return 0 cost 959ms``````

``````    Each sample counts as 0.01 seconds.
%   cumulative   self              self     total
time   seconds   seconds    calls  us/call  us/call  name
63.55      2.58     2.58   800000     3.23     3.23  is_prime(int)
15.76      3.22     0.64   800005     0.80     0.84  fermat_expmod(long, long, long)
10.10      3.63     0.41   800012     0.51     0.88  mr_expmod(long, long, long)
7.27      3.92     0.29 14813786     0.02     0.02  chk_unusual_square_root(long, long)
1.23      3.98     0.05                             main
0.86      4.01     0.04   800000     0.04     0.93  mr_test(long, long)
0.74      4.04     0.03 14813658     0.00     0.00  square(long)
0.25      4.05     0.01   800000     0.01     0.94  mr_is_prime(int)
0.25      4.06     0.01   800000     0.01     0.85  fermat_is_prime(int)
0.00      4.06     0.00   800000     0.00     0.84  fermat_test(long, long)``````