素数检测杂谈

最初我们学到的是这种朴素的写法：

(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 实现的，事先搞一个大数组，每一轮迭代都遍历数组，将已知素数的整数倍所对应的数组下标都划掉，最后留在数组中的就是素数。其实这个算法也是 O(N^2) 复杂度，随着 n 的增长，其代价也将迅速地变得不可接受。下面是我写的 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 上就有一个实现，顺便贡献了一个 O(logN) 的快速乘方取模算法：

(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 大概真的是个素数了。与费马检查不同，费马检查遇到 Carmichael 数无论重复多少次都无法给出正确的结果的。这就是它能避免 Carmichael 数愚弄的原因。

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

上面的代码我将迭代次数设为10次，事实上我用5次迭代来测试所有已知的 Carmichael 数，重复了很多次都没有发生错误。10次的迭代次数大概就可以把错误机率降到宇宙射线引发 CPU 执行指令错误的程度了吧。下面是测试代码和数表：

(define fname "Carmichael-list")

(call-with-input-file fname
  (lambda (ip)
    (let iter ((n (read ip)))
      (cond
       ((eof-object? n) #t)
       ((prime? n)
        (print n " test failed!" nl)
        (iter (read ip)))
       (else
        (iter (read ip)))))))

前方高能的分隔线

---

Carmichael-list:

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
530881
552721
656601
658801
670033
748657
825265
838201
852841
997633
1024651
1033669
1050985
1082809
1152271
1193221
1461241
1569457
1615681
1773289
1857241
1909001
2100901
2113921
2433601
2455921
2508013
2531845
2628073
2704801
3057601
3146221
3224065
3581761
3664585
3828001
4335241
4463641
4767841
4903921
4909177
5031181
5049001
5148001
5310721
5444489
5481451
5632705
5968873
6049681
6054985
6189121
6313681
6733693
6840001
6868261
7207201
7519441
7995169
8134561
8341201
8355841
8719309
8719921
8830801
8927101
9439201
9494101
9582145
9585541
9613297
9890881
10024561
10267951
10402561
10606681
10837321
10877581
11119105
11205601
11921001
11972017
12261061
12262321
12490201
12945745
13187665
13696033
13992265
14469841
14676481
14913991
15247621
15403285
15829633
15888313
16046641
16778881
17098369
17236801
17316001
17586361
17812081
18162001
18307381
18900973
19384289
19683001
20964961
21584305
22665505
23382529
25603201
26280073
26474581
26719701
26921089
26932081
27062101
27336673
27402481
28787185
29020321
29111881
31146661
31405501
31692805
32914441
33596641
34196401
34657141
34901461
35571601
35703361
36121345
36765901
37167361
37280881
37354465
37964809
38151361
38624041
38637361
39353665
40280065
40430401
40622401
40917241
41298985
41341321
41471521
42490801
43286881
43331401
43584481
43620409
44238481
45318561
45877861
45890209
46483633
47006785
48321001
49333201
50201089
53245921
54767881
55462177
56052361
58489201
60112885
60957361
62756641
64377991
64774081
65037817
65241793
67371265
67653433
67902031
67994641
68154001
69331969
70561921
72108421
72286501
74165065
75151441
75681541
75765313
76595761
77826001
78091201
78120001
79411201
79624621
80282161
80927821
81638401
81926461
82929001
83099521
83966401
84311569
84350561
84417985
87318001
88689601
90698401
92625121
93030145
93614521
93869665
94536001
96895441
99036001
99830641
99861985
100427041
101649241
101957401
102090781
104404861
104569501
104852881
105117481
105309289
105869401
107714881
109393201
109577161
111291181
114910489
115039081
115542505
116682721
118901521
119327041
120981601
121247281
122785741
124630273
127664461
128697361
129255841
129762001
130032865
130497361
132511681
133205761
133344793
133800661
134809921
134857801
135556345
136625941
139592101
139952671
140241361
144218341
145124785
146843929
150846961
151530401
151813201
153927961
157731841
158404141
158864833
159492061
161035057
161242705
161913961
163954561
167979421
168659569
169057801
169570801
170947105
171679561
172290241
172430401
172947529
173085121
174352641
175997185
176659201
178451857
178482151
178837201
180115489
181154701
182356993
184353001
186393481
186782401
188516329
188689501
189941761
193910977
194120389
194675041
196358977
200753281
206955841
208969201
212027401
214850881
214852609
216821881
221884001
226509361
227752993
228842209
230630401
230996949
231194965
237597361
238244041
238527745
241242001
242641153
246446929
247095361
250200721
252141121
255160621
256828321
257495641
258634741
266003101
270857521
271481329
271794601
273769921
274569601
275283401
277241401
278152381
279377281
280067761
288120421
289860481
291848401
292244833
292776121
295643089
295826581
296559361
299736181
300614161
301704985
302751505
306871201
311388337
321197185
321602401
328573477
329769721
333065305
333229141
338740417
334783585
346808881
348612265
354938221
357380101
358940737
360067201
362569201
364590721
366532321
366652201
367804801
367939585
368113411
382304161
382536001
390489121
392099401
393513121
393716701
395044651
395136505
399906001
403043257
405739681
413058601
413138881
416964241
419520241
426821473
429553345
434330401
434932961
438359041
440306461
455106601
458368201
461502097
461854261
462199681
471905281
471441001
473847121
477726145
481239361
483006889
484662529
490099681
490503601
492559141
503758801
507726901
510825601
511338241
516684961
517937581
518117041
518706721
527761081
529782121
530443201
532758241
540066241
542497201
544101481
545363281
547652161
548871961
549333121
549538081
551672221
552894301
555465601
556199281
556450777
557160241
558977761
561777121
564651361
568227241
569332177
573896881
577240273
579606301
580565233
590754385
595405201
597717121
600892993
602074585
602426161
606057985
609865201
612816751
616463809
620169409
625060801
625482001
629692801
631071001
633639097
652969351
656187001
662086041
683032801
683379841
686059921
689880801
697906561
702683101
703995733
704934361
710382401
710541481
711374401
713588401
717164449
727083001
739444021
743404663
744866305
752102401
759472561
765245881
771043201
775368901
775866001
776176261
784966297
790020001
790623289
794937601
798770161
804978721
809702401
809883361
814056001
822531841
824389441
829678141
833608321
834244501
839275921
841340521
843704401
847491361
849064321
851703301
851934601
852729121
855734401
863984881
867800701
876850801
882796321
885336481
888700681
897880321
902645857
914801665
918661501
928482241
931694401
934784929
935794081
939947009
940123801
941056273
954732853
955134181
957044881
958735681
958762729
958970545
962442001
962500561
963168193
968553181
975303121
977892241
981567505
981789337
985052881
990893569
993420289
993905641
1001152801
1027334881
1030401901
1031750401
1035608041
1038165961
1055384929
1070659201
1072570801
1093916341
1100674561
1103145121
1125038377
1131222841
1136739745
1177195201
1180398961
1189238401
1190790721
1193229577
1198650961
1200456577
1200778753
1207252621
1213619761
1216631521
1223475841
1227220801
1227280681
1232469001
1251295501
1251992281
1256855041
1257102001
1260332137
1264145401
1268604001
1269295201
1295577361
1299963601
1309440001
1312114945
1312332001
1316958721
1317828601
1318126321
1321983937
1332521065
1337805505
1348964401
1349671681
1376844481
1378483393
1382114881
1384157161
1394746081
1394942473
1404111241
1407548341
1422477001
1428966001
1439328001
1439492041
1441316269
1442761201
1490078305
1504651681
1507746241
1515785041
1520467201
1528936501
1540454761
1574601601
1576826161
1583582113
1588247851
1597821121
1626167341
1632785701
1646426881
1648076041
1659935761
1672719217
1676203201
1685266561
1688214529
1689411601
1690230241
1699279441
1701016801
1708549501
1726372441
1746692641
1750412161
1760460481
1772267281
1776450565
1778382541
1785507361
1795216501
1801558201
1803278401
1817067169
1825568641
1828377001
1831048561
1833328621
1841034961
1846817281
1848681121
1849811041
1879480513
1894344001
1899525601
1913016001
1918052065
1942608529
1943951041
1949646601
1950276565
1954174465
1955324449
1958102641
1976295241
1984089601
1988071801
2000436751
2023528501
2049293401
2064236401
2064373921
2067887557
2073560401
2080544005
2097317377
2101170097
2105594401
2107535221
2126689501
2140538401
2140699681
2301745249
2323147201
2436691321
2438403661
2444950561
2456536681
2529410281
2560600351
2598933481
2690867401
3264820001
3313196881
3480174001
3504570301
3713287801
3787491457
3835537861
4034969401
4215885697
4464573049
4562359201
4765950001
5255104513
5278692481
5781222721
5959748521
6047866621
6630702121
6916775113
7045248121
7629221377
8044161481
8251854001
8346731851
8652633601
8714965001
8976678481
9030158341
9086767201
9139810681
9293756581
9624742921
9701285761
9789244921
10119879001
10277275681
11346205609
12121569601
12173703001
12456671569
13079177569
13691148241
13946829751
14313548881
15906120889
16157879263
16765869121
17930792881
18151032901
18500666251
19742849041
21515221081
21796387201
22062297601
23224518901
23707055737
24285378001
25509692041
25749237001
26624905201
27278026129
30923424001
31876135201
34153717249
45983665729
48874811311
51436355851
51476019409
52756389001
57274147841
58094662081
59512667761
63593140801
64188507241
65700513721
67495942201
71171308081
82380774001
83946864769
100264053529
110296864801
168003672409
172018713961
173032371289
192739365541
225593397919
461574735553
464052305161
2199733160881
10028704049893
84154807001953
197531244744661
973694665856161
1436697831295441
1493812621027441
2094319836529921
2842648863161185
5778659093725441
6665161459969441
8015398984895401
8084842432668001
8188730132744161
12300849473799601
16166305446157945
32769125985828961
36629326277622001
39108343499765281
34795784213714161
37244219285276641
51116737346161201
61754693922936001
74538625278452401
91010482208190721
103183537035680641
126708084584398801
166857847951798081
186551892579535561
190232114399046721
194379756103868401
314610088213970641
335642734654849441
346413738355448401
556237362582392401
790689421836863641
1222628719906994401
1228155123631509217
1719756626091706801
2448237906710996401
2851896013395343201
3589102252820237401
4954039956700380001

转载于:https://www.cnblogs.com/zh-geek/p/6264774.html