素数检测杂谈
最初我们学到的是这种朴素的写法:
(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 上就有一个实现:
(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
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
· 探究高空视频全景AR技术的实现原理
· 理解Rust引用及其生命周期标识(上)
· DeepSeek 开源周回顾「GitHub 热点速览」
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布:重大改进与新特性概览!
· AI与.NET技术实操系列(二):开始使用ML.NET
· 单线程的Redis速度为什么快?