最开始我们学到的是这种朴素的写法:
(define (prime? n) (cond [(or (= n 2) (= n 3)) #t] [(even? n) #f] [else (prime-test-loop n)])) (define (prime-test-loop n) (let ((top (ceiling (sqrt n)))) (let iter ((start 3)) (cond [(> start top) #t] [(= (mod n start) 0) #f] [else (iter (+ start 2))]))))
虽然很笨拙,但其实这也是经过“优化”的,首先,除了 2 以外的偶数就不用判断了。其次,试除迭代从 3 起步,每次 +2,同样避开了偶数。最后,循环的结束条件是 (sqrt n).
后来又看到一种生成素数序列的方法,大体思想是维护一个已知的素数表,每一个新的数字 n 都用已知素数表里的数去试除。如果都不能整除,说明 n 是素数,将 n 添加到已知的素数表里,进入下一轮迭代。原文是用 C 实现的,事先搞一个大数组,每一轮迭代都遍历数组,将已知素数的整数倍所对应的数组下标都划掉,最后留在数组中的就是素数。下面是我写的 Scheme 版:
(define (prime-list n) (let iter ((start 3) (primes ‘(2))) (cond [(> start n) primes] [(andmap (lambda (p) (not (= (mod start p) 0))) primes) (iter (+ start 2) (cons start primes))] [else (iter (+ start 2) primes)])))
后来知道有一种基于概率的检测方法叫费马检查,SICP 上就有一个 Scheme 实现:
(define (square n) (* n n)) (define (expmod base exp m) (cond ((= exp 0) 1) ((even? exp) (mod (square (expmod base (/ exp 2) m)) m)) (else (mod (* base (expmod base (- exp 1) m)) m)))) (define (fermat-test n) (define (try-it a) (= (expmod a n n) a)) (try-it (+ 1 (random (- n 1))))) (define (fast-prime? n times) (cond ((= times 0) #t) ((fermat-test n) (fast-prime? n (- times 1))) (else #f)))
这个算法很不可靠,因为它所依据的理论是个伪命题,费马小定理断言:如果 n 是一个素数,那么小于 n 的任意正整数 a 的 n 次方再对 n 取模,结果依然为 a。
它的逆命题是,a 是任意一个小于 n 的正整数, 如果 a ^ n % n == a,那么 n 就是一个素数。很遗憾,这个逆命题不成立,因为有一类数叫做 Carmichael 数,它符合 a ^ n % n == a 这个特性,但却不是素数。直接用这段程序检查一个数是不是素数有很大的机会出错。这是已知 Carmichael 数的一部分:
561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461, 252601, 278545, 294409, 314821, 334153, 340561, 399001, 410041, 449065, 488881, 512461
可以看到,最小的数是561。如果你准备写一个程序,生成一个素数表,然而它迭代了500多次以后就开始胡言乱语,那玩笑就有点开大了。SICP 提到一个费马检测的变形,叫做 Miller-Rabin 检测,它能够避免 Carmichael 数的愚弄。
这个检测依据的断言是:如果 n 是一个素数,a 是任意一个小于 n 的正整数,那么 a ^ (n-1) % n == 1 % n。
要用 米勒-拉宾 检测算法测试一个数字是不是素数,我们应当选择一个小于 n 的随机正整数 a, 然后利用 expmod 函数计算 a ^ (n-1) % n . 但是当执行到 expmod 里的 square 那一步时,我们检测是不是找到了一个 (1 % n) 的非平凡的平方根,换言之,看一看是不是有这样一个数字,它既不是 1, 也不是 n-1,但是它的平方等于 1 % n。很容易证明,如果这样一个 1 的非平凡的平方根存在,那么 n 肯定不是素数。同样可以证明,如果 n 是一个并非素数的奇数,在小于 n 的正整数中,至少有一半在用这种方法计算 a ^ (n-1) 时能够找到这种非凡平方根。
中文版的 SICP 看不太明白,我自己翻译英文题目,还是觉得看不太明白。但大概意思明白了,如果一个小于 n 的正整数 a,它符合 a ^ (n-1) & n == 1,但是 a 既不是 1, 也不是 n-1,那么这个 n 肯定不是素数。对一个伪素数来说,只做一次测试,正确与错误的概率对半分,基本和算命差不多。如果做两次测试,那就很有可能发现它的真面目。重复的次数越多,结果越可信。如果重复多次依然通过了测试,那个这个数 n 大概真的是个素数了。这就是它能避免 Carmicheal 数愚弄的原因。
SICP 留了一个作业 1.28 ,要求实现 Miller-Rabin 检测,下面就是实现:
(define (expmod base exp m) (cond ((= exp 0) 1) ((even? exp) (check-nontrivial-sqrt (expmod base (/ exp 2) m) m)) ;; look here (else (mod (* base (expmod base (- exp 1) m)) m)))) (define (check-nontrivial-sqrt n m) (let ((x (mod (square n) m))) (if (and (not (= n 1)) (not (= n (- m 1))) (= x 1)) 0 x))) (define (miller-rabin-test n) (define (try-it a) (= (expmod a (- n 1) n) 1)) (try-it (+ 1 (random (- n 1))))) (define (fast-prime? n times) (cond ((= times 0) #t) ((miller-rabin-test n) (fast-prime? n (- times 1))) (else #f))) ;; 搞了个包装函数,偶数就不必判断了吧 (define (prime? n) (cond ((= n 2) #t) ((even? n) #f) (else (fast-prime? n 20))))
好吧,写到是写出来了,但是测试发现,它依然胡言乱语。。。
等有空再调试一下。