在 SICP 1,3 节的时候有提到判断素数的方法。
朴素的素数判定，地球的同学都知道的，如果数 n 除 1 以外的最小因数是他自己，那么n就是素数。
(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 )))
但这个方法很慢，它的阶是O(n)。
下面 SICP 给出了费马小定理：
如果一个数 n 是一个素数，而对而对任意 a<n，有 a^n≡a(mod n)
可惜哈，它的逆定理不成立，不过不代表不可用。我们可以挑选小于 n 的数 a，做上文所述的测试。如果满足 a^n≡a(mod n)，那么 n 很可能是一个素数。而挑选越多的 a 对 n 进行测试，n 是素数的可能性也越大。所以这是一个人品算法，只要测试的数量足够，那么 n 不是素数的概率就会变得很小。
于是我们有了 fermat 测试的程序。不用担心求 a^n 的效率问题，因为比较好解决：
试看：
a^n = (a^(n/2))^2; 当n为偶数
a^n = a* a^(n-1); 当n为奇数
总是成立的，所以几乎每次作乘法的时候，n 就增加了一倍。放到数轴看一看，其实这是二分法。所以时间复杂度得到降低。所以大概知道为什么编程语言提供了求 n 次幂的函数 pow 总是比自己写的 a*a*a… 快了。
但是不应该直接使用语言自带的 pow 函数。原因是我们的目标不是 a 的 n 次幂，而是想知道 a 的 n 次幂对 m 取模是否等于 a 对 m 取模。简化一下，即只需判断这个表达式是否成立：a^n mod m = a。如果直接表达成 pow(a, n) % m == a 计算机会郁闷，因为 pow(a, n) 的值会很大。所以应该采用类似的技术：
a^n mod m = ((a^(n / 2))^2) mod m; 当n为偶数
a^n mod m = (a* a^(n-1)) mod m; 当n为奇数
这样，每次计算得到的值都不比 m 大了，避免过大的数参与计算。于是很容易有了以下程序
(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)))
但是，总有一些不完美。
存在一种数叫 Carmichael 数，能骗过费马测试。在 2 到 10000 这个区间内有这么几个：
561，1105，1729，2465，2821，6601，8911
这些钉子户是合数，但是不管你用什么 a 去测试，他们都可以在费马测试中冒充素数。他们数量也不少，而分布广，在100 000 000 以内一共有 255 个 Carmichael 数。所以，如果不愿意去记录 Carmichael 数表的话，可以加强我们的素数测试方法。
还好，在练习 1.28中，Miller 和 Rabin 给我们带来了 Miller-Rabin 测试法。该测试法基于 fermat 小定理的一个变形：
如果 n 是素数，对任何 a < n，a^(n-1) ≡ 1(mod n)
要使用 MR 测试法呢，就需要在执行上文 Fermat 测试的平方步骤时检查是否遇到了一个 “1 mod n的非平凡平方根”，如果存在一个 s = a^2; 而且 s <> 1 且 s <> n-1，而且 s^2 mod m = 1，那么这个数n就不是素数了。
Matrix67：Miller-Rabin 素性测试同样是不确定算法，我们把可以通过以a为底的 Miller-Rabin 测试的合数称作以a为底的强伪素数(strong pseudoprime)。第一个以 2 为底的强伪素数为 2047。第一个以 2 和 3 为底的强伪素数则大到 1 373 653。
所以当你测试量足够，那么该方法可以在工业应用中保证安全，下面是我写的1.28题答案
(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)))
我构造了一些数据，分别求[1000, 1100], [10000, 10100], [100000, 100100], [1000000, 1000100]这些区间的素数，看看它们这些方法的表现如何：
常规方法 常规方法的因数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
这是在我的 intel Core 1 1.83GHz 的机器上跑的结果（Ubuntu 8.04 + DrScheme）。
我还另写了一个 C 语言的版本，运行结果如下：
区间[2, 10000]
$ ./prime
Normal: return 0 cost 10ms
Fermat: return 0 cost 24ms
Miller-Rabin: return 0 cost 28ms
区间[200000, 1000000]:
$ ./prime
Normal: return 0 cost 3088ms
Fermat: return 0 cost 824ms
Miller-Rabin: return 0 cost 959ms
采用 gprof 分析的结果基本与预期一致，部分分析报告如下：
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)
当 n 比较小的时候效果不好的原因我猜是 C 语言在处理普通递归时花销比较大。而当数据变得庞大时，这种递归上的花销的影响就相对不大了。如果有同学把上面的算法改成迭代形式（我不会改:P）或者是尾递归形式（gcc可以把尾递归形式优化成迭代形式）我想效果更加明显。