SICP虽然已经拿到手很长时间了，一直是懒洋洋地看，快看到完第3章了，就是没做题。昨天遇到了有关判定素数的问题，似乎见过的，决定把前面的习题做了….我错了，:( 然后我决定以后的习题每题都要做。

好吧。在1,3节的时候有提到判断素数的方法。朴素的素数判定，地球的同学都知道的，如果数n除1以外的最小因数是他自己，那么n就是素数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14

 
;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)的，当n变得很大的时候，会很慢。
下面SICP给出了费马小定理：

如果一个数n是一个素数，而对任意a<n，有a^n≡a(mod n)

可惜哈，它的逆定理是不成立di，不过不代表不可用。我们可以挑选小于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大了，避免过分大的数参与计算。于是很容易有了一下程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

 
(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)))

但是（我们一般最怕MM在关键时刻说“但是”），总有一些不完美。那就是存在一种数长得欠扁，叫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。

所以当你测试量足够，那么该方法可以在工业应用保证安全（看上去SICP个理论派，其实还是实践派的，NND，在MIT计算机属于工学，我还在看他们的入门教材，NND），下面是我写的1.28题答案

 
(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分析的结果基本上和预期一致，这是部分分析报告：

Flat profile:

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可以把尾递归形式优化成迭代形式）我想效果更加明显。

Related posts(相关文章)

• 对付冗余计算…的那个途径 (comment: 14)
• Church计数和Lambda演算 (comment: 2)
• 阿西莫夫机器人定律与软件设计 (comment: 4)
• 这篇日志是我在微软中国比赛时发表的 (comment: 2)
• 自己做一个cairo clock (comment: 1)