[原创]浅析求素数算法

chai2010

32位无符号乘法，返回高32位。
低32位可以直接(u*v)得到

1. 
   
2. unsigned mulhu(unsigned u, unsigned v)
   
3. {
   
4.    unsigned u0, v0, w0;
   
5.    unsigned u1, v1, w1, w2, t;
   
6. 
   
7.    u0 = u & 0xFFFF;  u1 = u >> 16;
   
8.    v0 = v & 0xFFFF;  v1 = v >> 16;
   
9.    w0 = u0*v0;
   
10.    t  = u1*v0 + (w0 >> 16);
    
11.    w1 = t & 0xFFFF;
    
12.    w2 = t >> 16;
    
13.    w1 = u0*v1 + w1;
    
14. 
    
15.    return u1*v1 + w2 + (w1 >> 16);
    
16. }
    

复制代码

chai2010

(x*y)%m代码修改如下

1. 
   
2. 
   
3. // 计算(u*v)%m
   
4. 
   
5. unsigned mul_mod(unsigned u, unsigned v, unsigned z)
   
6. {
   
7.         // 如果u*v没有溢出, 则直接计算
   
8. 
   
9.         if(u*v >= 1.0*u*v) return (u*v)%z;
   
10. 
    
11.         // 进行长乘法(结果为64位)
    
12. 
    
13.         unsigned u0, v0, w0;
    
14.         unsigned u1, v1, w1, w2, t;
    
15. 
    
16.         u0 = u & 0xFFFF;  u1 = u >> 16;
    
17.         v0 = v & 0xFFFF;  v1 = v >> 16;
    
18.         w0 = u0*v0;
    
19.         t  = u1*v0 + (w0 >> 16);
    
20.         w1 = t & 0xFFFF;
    
21.         w2 = t >> 16;
    
22.         w1 = u0*v1 + w1;
    
23. 
    
24.         // x为高32位, y为低32位
    
25. 
    
26.         unsigned x = u1*v1 + w2 + (w1 >> 16);
    
27.         unsigned y = u*v;
    
28. 
    
29.         // 进行长除法(除数为64位)
    
30. 
    
31.         for (int i = 1; i <= 32; i++)
    
32.         {
    
33.                 t = (int)x >> 31;                // All 1's if x(31) = 1.
    
34. 
    
35.                 x = (x << 1) | (y >> 31);        // Shift x || y left
    
36.                 y <<= 1;                        // one bit.
    
37.                
    
38.                 if((x|t) >= z) { x -= z; y++; }
    
39.         }
    
40. 
    
41.         return x;        // y为商, x为余数
    
42. }
    
43. 
    
44. // 判断是否是费马(伪)素数
    
45. 
    
46. static bool isFermatNum(unsigned p)
    
47. {
    
48.         unsigned m = p--;
    
49.         unsigned r = 2%m;
    
50.         unsigned k = 1;
    
51.    
    
52.         // 蒙格马利快速幂模算法
    
53. 
    
54.         while(p > 1)
    
55.         {
    
56.                 if(p&1) k = mul_mod(k, r, m);
    
57.                 r = mul_mod(r, r, m); p >>= 1;
    
58.         }
    
59. 
    
60.         return mul_mod(k, r, m) == 1;
    
61. }
    

复制代码

chai2010

整理之后

1. 
   
2. 
   
3. // 函数: bool isPrime(unsigned n);
   
4. // 描述: 判断n(0, 2^32-1)是否为素数
   
5. // 算法: 费马小定理-卡尔麦克数
   
6. // 作者: 柴树杉
   
7. // 时间: 2006-11-06
   
8. // 主页: [url]chaishushan.googlepages.com[/url]
   
9. 
   
10. #include <assert.h>
    
11. #include <stdlib.h>
    
12. 
    
13. // 计算数组的元素个数
    
14. 
    
15. #define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))
    
16. 
    
17. // 2^32范围内的所有卡尔麦克数
    
18. 
    
19. static const unsigned CarmichaeNumbers[] =
    
20. {
    
21. 561,1105,1729,2465,2821,6601,8911,10585,15841,29341,41041,46657,52633,62745,
    
22. 63973,75361,101101,115921,126217,162401,172081,188461,252601,278545,294409,
    
23. 314821,334153,340561,399001,410041,449065,488881,512461,530881,552721,656601,
    
24. 658801,670033,748657,825265,838201,852841,997633,1024651,1033669,1050985,
    
25. 1082809,1152271,1193221,1461241,1569457,1615681,1773289,1857241,1909001,
    
26. 2100901,2113921,2433601,2455921,2508013,2531845,2628073,2704801,3057601,
    
27. 3146221,3224065,3581761,3664585,3828001,4335241,4463641,4767841,4903921,
    
28. 4909177,5031181,5049001,5148001,5310721,5444489,5481451,5632705,5968873,
    
29. 6049681,6054985,6189121,6313681,6733693,6840001,6868261,7207201,7519441,
    
30. 7995169,8134561,8341201,8355841,8719309,8719921,8830801,8927101,9439201,
    
31. 9494101,9582145,9585541,9613297,9890881,10024561,10267951,10402561,10606681,
    
32. 10837321,10877581,11119105,11205601,11921001,11972017,12261061,12262321,
    
33. 12490201,12945745,13187665,13696033,13992265,14469841,14676481,14913991,
    
34. 15247621,15403285,15829633,15888313,16046641,16778881,17098369,17236801,
    
35. 17316001,17586361,17812081,18162001,18307381,18900973,19384289,19683001,
    
36. 20964961,21584305,22665505,23382529,25603201,26280073,26474581,26719701,
    
37. 26921089,26932081,27062101,27336673,27402481,28787185,29020321,29111881,
    
38. 31146661,31405501,31692805,32914441,33302401,33596641,34196401,34657141,
    
39. 34901461,35571601,35703361,36121345,36765901,37167361,37280881,37354465,
    
40. 37964809,38151361,38624041,38637361,39353665,40160737,40280065,40430401,
    
41. 40622401,40917241,41298985,41341321,41471521,42490801,43286881,43331401,
    
42. 43584481,43620409,44238481,45318561,45877861,45890209,46483633,47006785,
    
43. 48321001,48628801,49333201,50201089,53245921,53711113,54767881,55462177,
    
44. 56052361,58489201,60112885,60957361,62756641,64377991,64774081,65037817,
    
45. 65241793,67371265,67653433,67902031,67994641,68154001,69331969,70561921,
    
46. 72108421,72286501,74165065,75151441,75681541,75765313,76595761,77826001,
    
47. 78091201,78120001,79411201,79624621,80282161,80927821,81638401,81926461,
    
48. 82929001,83099521,83966401,84311569,84350561,84417985,87318001,88689601,
    
49. 90698401,92625121,93030145,93614521,93869665,94536001,96895441,99036001,
    
50. 99830641,99861985,100427041,101649241,101957401,102090781,104404861,104569501,
    
51. 104852881,105117481,105309289,105869401,106041937,107714881,109393201,109577161,
    
52. 111291181,114910489,115039081,115542505,116682721,118901521,119327041,120981601,
    
53. 121247281,122785741,124630273,127664461,128697361,129255841,129762001,130032865,
    
54. 130497361,132511681,133205761,133344793,133800661,134809921,134857801,135556345,
    
55. 136625941,139592101,139952671,140241361,144218341,145124785,146843929,150846961,
    
56. 151530401,151813201,153589801,153927961,157731841,158404141,158864833,159492061,
    
57. 161035057,161242705,161913961,163954561,167979421,168659569,169057801,169570801,
    
58. 170947105,171454321,171679561,172290241,172430401,172947529,173085121,174352641,
    
59. 175997185,176659201,178451857,178482151,178837201,180115489,181154701,182356993,
    
60. 184353001,186393481,186782401,187188001,188516329,188689501,189941761,193708801,
    
61. 193910977,194120389,194675041,196358977,200753281,206955841,208969201,212027401,
    
62. 213835861,214850881,214852609,216821881,221884001,225745345,226509361,227752993,
    
63. 228842209,230630401,230996949,231194965,237597361,238244041,238527745,241242001,
    
64. 242641153,246446929,247095361,250200721,252141121,255160621,256828321,257495641,
    
65. 258634741,266003101,270857521,271481329,271794601,273769921,274569601,275283401,
    
66. 277241401,278152381,279377281,280067761,280761481,288120421,289766701,289860481,
    
67. 291848401,292244833,292776121,295643089,295826581,296559361,299736181,300614161,
    
68. 301704985,302751505,306871201,311388337,318266641,321197185,321602401,325546585,
    
69. 328573477,329769721,333065305,333229141,334783585,338740417,346808881,348612265,
    
70. 354938221,357277921,357380101,358940737,360067201,362569201,364590721,366532321,
    
71. 366652201,367804801,367939585,368113411,382304161,382536001,390489121,392099401,
    
72. 393122521,393513121,393716701,395044651,395136505,396262945,399906001,403043257,
    
73. 403317421,405739681,413058601,413138881,413631505,416964241,417241045,419520241,
    
74. 426821473,429553345,434330401,434932961,438359041,440306461,440707345,455106601,
    
75. 458368201,461502097,461854261,462199681,471441001,471905281,473847121,477726145,
    
76. 481239361,483006889,484662529,490099681,490503601,492559141,496050841,499310197,
    
77. 503758801,507726901,509033161,510825601,511338241,516684961,517937581,518117041,
    
78. 518706721,527761081,529782121,530443201,532758241,533860309,540066241,542497201,
    
79. 544101481,545363281,545570641,547652161,548871961,549117205,549333121,549538081,
    
80. 551672221,552894301,555465601,556199281,556450777,557160241,557795161,558570961,
    
81. 558977761,561481921,561777121,564651361,568227241,569332177,573896881,577240273,
    
82. 579606301,580565233,590754385,593234929,595405201,597717121,600892993,602074585,
    
83. 602426161,602593441,606057985,609865201,611397865,612347905,612816751,616463809,
    
84. 618068881,620169409,621101185,625060801,625482001,629692801,631071001,633639097,
    
85. 638959321,642708001,652969351,656187001,662086041,672389641,683032801,683379841,
    
86. 684106401,686059921,689880801,697906561,698548201,702683101,703995733,704934361,
    
87. 705101761,707926801,710382401,710541481,711374401,713588401,717164449,721244161,
    
88. 722923201,727083001,739444021,743404663,744866305,745864945,746706961,752102401,
    
89. 759472561,765245881,771043201,775368901,775866001,776176261,781347841,784966297,
    
90. 790020001,790623289,794937601,798770161,804978721,809702401,809883361,814056001,
    
91. 822531841,824389441,824405041,829678141,832060801,833608321,834244501,834720601,
    
92. 836515681,839275921,841340521,842202361,843704401,847491361,849064321,851703301,
    
93. 851934601,852729121,854197345,855734401,860056705,863984881,867800701,868234081,
    
94. 876850801,882796321,885336481,888700681,891706861,897880321,902645857,914801665,
    
95. 918661501,928482241,930745621,931694401,934784929,935794081,939947009,940123801,
    
96. 941056273,945959365,947993761,954732853,955134181,957044881,958735681,958762729,
    
97. 958970545,962442001,962500561,963163201,963168193,968553181,968915521,975303121,
    
98. 977737321,977892241,981567505,981789337,985052881,986088961,990893569,993420289,
    
99. 993905641,1001152801,1018928485,1027334881,1030401901,1031750401,1035608041,
    
100. 1038165961,1055384929,1070659201,1072570801,1074363265,1079556193,1090842145,
     
101. 1093916341,1100674561,1103145121,1125038377,1131222841,1132988545,1134044821,
     
102. 1136739745,1138049137,1140441121,1150270849,1152793621,1162202581,1163659861,
     
103. 1177195201,1177800481,1180398961,1183104001,1189238401,1190790721,1193229577,
     
104. 1194866101,1198650961,1200456577,1200778753,1206057601,1207252621,1210178305,
     
105. 1213619761,1214703721,1216631521,1223475841,1227220801,1227280681,1232469001,
     
106. 1251295501,1251992281,1254318481,1256855041,1257102001,1260332137,1264145401,
     
107. 1268604001,1269295201,1271325841,1295577361,1299963601,1309440001,1312114945,
     
108. 1312332001,1316958721,1317828601,1318126321,1321983937,1330655041,1332521065,
     
109. 1337805505,1348964401,1349671681,1362132541,1376844481,1378483393,1382114881,
     
110. 1384157161,1394746081,1394942473,1404111241,1404228421,1407548341,1410833281,
     
111. 1419339691,1420379065,1422477001,1423668961,1428966001,1439328001,1439492041,
     
112. 1441316269,1442761201,1445084173,1452767521,1481619601,1483199641,1490078305,
     
113. 1490522545,1504651681,1505432881,1507746241,1515785041,1520467201,1521835381,
     
114. 1528936501,1529544961,1534274841,1540454761,1545387481,1571503801,1573132561,
     
115. 1574601601,1576826161,1583582113,1588247851,1592668441,1597009393,1597821121,
     
116. 1615335085,1618206745,1626167341,1632785701,1641323905,1646426881,1648076041,
     
117. 1657700353,1659935761,1672719217,1674309385,1676203201,1676641681,1678569121,
     
118. 1685266561,1688214529,1688639041,1689411601,1690230241,1696572001,1698623641,
     
119. 1699279441,1701016801,1705470481,1708549501,1726372441,1746692641,1750412161,
     
120. 1752710401,1760460481,1769031901,1772267281,1773486001,1776450565,1778382541,
     
121. 1784291041,1784975941,1785507361,1795216501,1801558201,1803278401,1817067169,
     
122. 1825568641,1828377001,1831048561,1833328621,1836304561,1841034961,1845871105,
     
123. 1846817281,1848112761,1848681121,1849811041,1854001513,1855100017,1858395529,
     
124. 1879480513,1887933601,1894344001,1896961801,1899525601,1913016001,1916729101,
     
125. 1918052065,1919767681,1932608161,1942608529,1943951041,1949646601,1950276565,
     
126. 1954174465,1955324449,1958102641,1962804565,1976295241,1984089601,1988071801,
     
127. 1992841201,1999004365,2000436751,2023528501,2029554241,2048443501,2048751901,
     
128. 2049293401,2064236401,2064373921,2067887557,2073560401,2080544005,2097317377,
     
129. 2101170097,2105594401,2107535221,2111416021,2111488561,2117725921,2126689501,
     
130. 2140538401,2140699681,2159003281,2170282969,2176838049,2178944461,2199700321,
     
131. 2199931651,2201169601,2212935985,2215407601,2216430721,2217951073,2223876601,
     
132. 2224519921,2230305949,2232385345,2239622113,2244932281,2246916001,2258118721,
     
133. 2265650401,2272748401,2278677961,2295209281,2301745249,2302419601,2309027281,
     
134. 2313774001,2320224481,2320690177,2322648901,2323147201,2332627249,2335640077,
     
135. 2339165521,2342644921,2353639681,2359686241,2361232477,2367379201,2377166401,
     
136. 2391137281,2396357041,2407376665,2414829781,2430556381,2436691321,2438403661,
     
137. 2443829641,2444950561,2456536681,2457411265,2467813621,2470348441,2470894273,
     
138. 2479305985,2480343553,2489462641,2494465921,2494621585,2497638781,2509860961,
     
139. 2510085721,2519819281,2523947041,2527812001,2529410281,2539024741,2544590161,
     
140. 2547621973,2560104001,2560600351,2561945401,2564889601,2573686441,2574243721,
     
141. 2575260241,2586927553,2588653081,2597928961,2598933481,2601144001,2602378721,
     
142. 2605557781,2607162961,2607237361,2616662881,2617181281,2630374741,2642025673,
     
143. 2657502001,2665141921,2677147201,2685422593,2690867401,2693939401,2702470861,
     
144. 2709611521,2723859001,2733494401,2735309521,2766172501,2766901501,2770560241,
     
145. 2776874941,2787998641,2797002901,2801124001,2806205689,2811315361,2815304401,
     
146. 2832480001,2833846561,2842912381,2858298301,2867755969,2942952481,2943556201,
     
147. 2965085641,2998467901,3001561441,3007991701,3022354401,3024774901,3025708561,
     
148. 3030758401,3034203361,3035837161,3044238121,3044970001,3068534701,3069196417,
     
149. 3072080089,3072094201,3077802001,3078386641,3086434561,3088134721,3090578401,
     
150. 3102234751,3104207821,3105567361,3112974481,3119101921,3138302401,3159422785,
     
151. 3159939601,3164207761,3180288385,3180632833,3188744065,3190894201,3193414093,
     
152. 3203895601,3215031751,3222053185,3232450585,3240392401,3245477761,3246238801,
     
153. 3248891101,3249390145,3263564305,3264820001,3270933121,3277595665,3281736601,
     
154. 3284630713,3307322305,3313196881,3313744561,3314111761,3319323601,3328437481,
     
155. 3345878017,3347570941,3348463105,3353809537,3378014641,3380740301,3411338491,
     
156. 3413656441,3429457921,3438721441,3441837421,3480174001,3504570301,3508507801,
     
157. 3521441665,3534510001,3555636481,3574014445,3575798785,3576804001,3600918181,
     
158. 3618244081,3630291841,3637405045,3682471321,3697952401,3711456001,3712280041,
     
159. 3713287801,3715938721,3722793481,3727589761,3754483201,3767865601,3776698801,
     
160. 3787491457,3799111681,3800513761,3801823441,3805181281,3832646221,3834444901,
     
161. 3835537861,3858853681,3863326897,3880251649,3891338101,3892568065,3901730401,
     
162. 3907357441,3922752121,3928256641,3951813601,3981047941,3998554561,4015029061,
     
163. 4030864201,4034969401,4059151489,4065133501,4077957961,4115677501,4127050621,
     
164. 4134273793,4138747921,4146685921,4160472121,4162880401,4167038161,4169092201,
     
165. 4169867689,4189909501,4199202001,4199529601,4199932801,4202009461,4210922233,
     
166. 4212413569,4215885697,4216799521,4277982241        // 1119: 4295605861, 溢出
     
167. };
     
168. 
     
169. // hash表的改进方法, 只需要保存下标
     
170. // 经测试, 最多有6次冲突, b[i][0]保存数目
     
171. // 注意buckets的各维大小是人工测试得到!!
     
172. 
     
173. static short buckets[1021][6+1];
     
174. 
     
175. // 自动初始化hash表
     
176. 
     
177. static struct init
     
178. {
     
179.         init(void)
     
180.         {
     
181.                 for(int i = 0; i < NELEMS(CarmichaeNumbers); ++i)
     
182.                 {
     
183.                         // buckets[v]类似一个栈
     
184.                         // buckets[v][0]记录栈顶位置
     
185.                         // 最多有6次冲突, 不会溢出
     
186. 
     
187.                         int key = CarmichaeNumbers[i]%NELEMS(buckets);
     
188.                         assert(buckets[key][0] < NELEMS(buckets[0])-1);
     
189. 
     
190.                         buckets[key][++buckets[key][0]] = i;
     
191.                 }
     
192.         }
     
193. } _init;
     
194. 
     
195. // 判断是否是卡尔麦克数
     
196. 
     
197. static bool isCarmichaelNum(unsigned n)
     
198. {
     
199.         // 改进: 采用hash表查找
     
200. 
     
201.         int i, key = n%NELEMS(buckets);
     
202. 
     
203.         for(i = 1; i <= buckets[key][0]; ++i)
     
204.                 if(n == CarmichaeNumbers[buckets[key][i]]) return true;
     
205. 
     
206.         return false;
     
207. }
     
208. 
     
209. // 计算(u*v)%m
     
210. 
     
211. unsigned mul_mod(unsigned u, unsigned v, unsigned z)
     
212. {
     
213.         // 如果u*v没有溢出, 则直接计算
     
214. 
     
215.         if((u*v)/u == v) return (u*v)%z;
     
216. 
     
217.         // 进行长乘法(结果为64位)
     
218. 
     
219.         unsigned u0, v0, w0;
     
220.         unsigned u1, v1, w1, w2, t;
     
221. 
     
222.         u0 = u & 0xFFFF;  u1 = u >> 16;
     
223.         v0 = v & 0xFFFF;  v1 = v >> 16;
     
224.         w0 = u0*v0;
     
225.         t  = u1*v0 + (w0 >> 16);
     
226.         w1 = t & 0xFFFF;
     
227.         w2 = t >> 16;
     
228.         w1 = u0*v1 + w1;
     
229. 
     
230.         // x为高32位, y为低32位
     
231. 
     
232.         unsigned x = u1*v1 + w2 + (w1 >> 16);
     
233.         unsigned y = u*v;
     
234. 
     
235.         // 进行长除法(被除数64位, 除数32位)
     
236. 
     
237.         for (int i = 1; i <= 32; i++)
     
238.         {
     
239.                 t = (int)x >> 31;           // All 1's if x(31) = 1.
     
240. 
     
241.                 x = (x << 1) | (y >> 31);   // Shift x || y left
     
242.                 y <<= 1;                    // one bit.
     
243.                
     
244.                 if((x|t) >= z) { x -= z; y++; }
     
245.         }
     
246. 
     
247.         return x;        // y为商, x为余数
     
248. }
     
249. 
     
250. // 判断是否是费马(伪)素数
     
251. 
     
252. static bool isFermatNum(unsigned p)
     
253. {
     
254.         unsigned m = p--;
     
255.         unsigned r = 2%m;
     
256.         unsigned k = 1;
     
257.    
     
258.         // 蒙格马利快速幂模算法
     
259. 
     
260.         while(p > 1)
     
261.         {
     
262.                 if(p&1) k = mul_mod(k, r, m);
     
263.                 r = mul_mod(r, r, m); p >>= 1;
     
264.         }
     
265. 
     
266.         return mul_mod(k, r, m) == 1;
     
267. }
     
268. 
     
269. // 根据费马小定理测试
     
270. 
     
271. bool isPrime(unsigned n)
     
272. {
     
273.         // 处理比较特殊的情况
     
274. 
     
275.         if(n < 2) return false;
     
276.         if(n == 2) return true;
     
277.         if(!(n&1)) return false;
     
278. 
     
279.         // 是费马(伪)素数但不是卡尔麦克数则必定为素数
     
280. 
     
281.         return isFermatNum(n) && !isCarmichaelNum(n);
     
282. }
     

复制代码

chai2010

Post by chai2010
整理之后

1. 
   
2. 
   
3. // 函数: bool isPrime(unsigned n);
   
4. // 描述: 判断n(0, 2^32-1)是否为素数
   
5. // 算法: 费马小定理-卡尔麦克数
   
6. // 作者: 柴树杉
   
7. // 时间: 2006-11-06
   
8. // 主页: [url]chaishushan.googlepages.com[/url]
   
9. 
   
10. #include <assert.h>
    
11. #include <stdlib.h>
    
12. 
    
13. // 计算数组的元素个数
    
14. 
    
15. #define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))
    
16. 
    
17. // 2^32范围内的所有卡尔麦克数
    
18. 
    
19. static const unsigned CarmichaeNumbers[] =
    
20. {
    
21. 561,1105,1729,2465,2821,6601,8911,10585,15841,29341,41041,46657,52633,62745,
    
22. 63973,75361,101101,115921,126217,162401,172081,188461,252601,278545,294409,
    
23. 314821,334153,340561,399001,410041,449065,488881,512461,530881,552721,656601,
    
24. 658801,670033,748657,825265,838201,852841,997633,1024651,1033669,1050985,
    
25. 1082809,1152271,1193221,1461241,1569457,1615681,1773289,1857241,1909001,
    
26. 2100901,2113921,2433601,2455921,2508013,2531845,2628073,2704801,3057601,
    
27. 3146221,3224065,3581761,3664585,3828001,4335241,4463641,4767841,4903921,
    
28. 4909177,5031181,5049001,5148001,5310721,5444489,5481451,5632705,5968873,
    
29. 6049681,6054985,6189121,6313681,6733693,6840001,6868261,7207201,7519441,
    
30. 7995169,8134561,8341201,8355841,8719309,8719921,8830801,8927101,9439201,
    
31. 9494101,9582145,9585541,9613297,9890881,10024561,10267951,10402561,10606681,
    
32. 10837321,10877581,11119105,11205601,11921001,11972017,12261061,12262321,
    
33. 12490201,12945745,13187665,13696033,13992265,14469841,14676481,14913991,
    
34. 15247621,15403285,15829633,15888313,16046641,16778881,17098369,17236801,
    
35. 17316001,17586361,17812081,18162001,18307381,18900973,19384289,19683001,
    
36. 20964961,21584305,22665505,23382529,25603201,26280073,26474581,26719701,
    
37. 26921089,26932081,27062101,27336673,27402481,28787185,29020321,29111881,
    
38. 31146661,31405501,31692805,32914441,33302401,33596641,34196401,34657141,
    
39. 34901461,35571601,35703361,36121345,36765901,37167361,37280881,37354465,
    
40. 37964809,38151361,38624041,38637361,39353665,40160737,40280065,40430401,
    
41. 40622401,40917241,41298985,41341321,41471521,42490801,43286881,43331401,
    
42. 43584481,43620409,44238481,45318561,45877861,45890209,46483633,47006785,
    
43. 48321001,48628801,49333201,50201089,53245921,53711113,54767881,55462177,
    
44. 56052361,58489201,60112885,60957361,62756641,64377991,64774081,65037817,
    
45. 65241793,67371265,67653433,67902031,67994641,68154001,69331969,70561921,
    
46. 72108421,72286501,74165065,75151441,75681541,75765313,76595761,77826001,
    
47. 78091201,78120001,79411201,79624621,80282161,80927821,81638401,81926461,
    
48. 82929001,83099521,83966401,84311569,84350561,84417985,87318001,88689601,
    
49. 90698401,92625121,93030145,93614521,93869665,94536001,96895441,99036001,
    
50. 99830641,99861985,100427041,101649241,101957401,102090781,104404861,104569501,
    
51. 104852881,105117481,105309289,105869401,106041937,107714881,109393201,109577161,
    
52. 111291181,114910489,115039081,115542505,116682721,118901521,119327041,120981601,
    
53. 121247281,122785741,124630273,127664461,128697361,129255841,129762001,130032865,
    
54. 130497361,132511681,133205761,133344793,133800661,134809921,134857801,135556345,
    
55. 136625941,139592101,139952671,140241361,144218341,145124785,146843929,150846961,
    
56. 151530401,151813201,153589801,153927961,157731841,158404141,158864833,159492061,
    
57. 161035057,161242705,161913961,163954561,167979421,168659569,169057801,169570801,
    
58. 170947105,171454321,171679561,172290241,172430401,172947529,173085121,174352641,
    
59. 175997185,176659201,178451857,178482151,178837201,180115489,181154701,182356993,
    
60. 184353001,186393481,186782401,187188001,188516329,188689501,189941761,193708801,
    
61. 193910977,194120389,194675041,196358977,200753281,206955841,208969201,212027401,
    
62. 213835861,214850881,214852609,216821881,221884001,225745345,226509361,227752993,
    
63. 228842209,230630401,230996949,231194965,237597361,238244041,238527745,241242001,
    
64. 242641153,246446929,247095361,250200721,252141121,255160621,256828321,257495641,
    
65. 258634741,266003101,270857521,271481329,271794601,273769921,274569601,275283401,
    
66. 277241401,278152381,279377281,280067761,280761481,288120421,289766701,289860481,
    
67. 291848401,292244833,292776121,295643089,295826581,296559361,299736181,300614161,
    
68. 301704985,302751505,306871201,311388337,318266641,321197185,321602401,325546585,
    
69. 328573477,329769721,333065305,333229141,334783585,338740417,346808881,348612265,
    
70. 354938221,357277921,357380101,358940737,360067201,362569201,364590721,366532321,
    
71. 366652201,367804801,367939585,368113411,382304161,382536001,390489121,392099401,
    
72. 393122521,393513121,393716701,395044651,395136505,396262945,399906001,403043257,
    
73. 403317421,405739681,413058601,413138881,413631505,416964241,417241045,419520241,
    
74. 426821473,429553345,434330401,434932961,438359041,440306461,440707345,455106601,
    
75. 458368201,461502097,461854261,462199681,471441001,471905281,473847121,477726145,
    
76. 481239361,483006889,484662529,490099681,490503601,492559141,496050841,499310197,
    
77. 503758801,507726901,509033161,510825601,511338241,516684961,517937581,518117041,
    
78. 518706721,527761081,529782121,530443201,532758241,533860309,540066241,542497201,
    
79. 544101481,545363281,545570641,547652161,548871961,549117205,549333121,549538081,
    
80. 551672221,552894301,555465601,556199281,556450777,557160241,557795161,558570961,
    
81. 558977761,561481921,561777121,564651361,568227241,569332177,573896881,577240273,
    
82. 579606301,580565233,590754385,593234929,595405201,597717121,600892993,602074585,
    
83. 602426161,602593441,606057985,609865201,611397865,612347905,612816751,616463809,
    
84. 618068881,620169409,621101185,625060801,625482001,629692801,631071001,633639097,
    
85. 638959321,642708001,652969351,656187001,662086041,672389641,683032801,683379841,
    
86. 684106401,686059921,689880801,697906561,698548201,702683101,703995733,704934361,
    
87. 705101761,707926801,710382401,710541481,711374401,713588401,717164449,721244161,
    
88. 722923201,727083001,739444021,743404663,744866305,745864945,746706961,752102401,
    
89. 759472561,765245881,771043201,775368901,775866001,776176261,781347841,784966297,
    
90. 790020001,790623289,794937601,798770161,804978721,809702401,809883361,814056001,
    
91. 822531841,824389441,824405041,829678141,832060801,833608321,834244501,834720601,
    
92. 836515681,839275921,841340521,842202361,843704401,847491361,849064321,851703301,
    
93. 851934601,852729121,854197345,855734401,860056705,863984881,867800701,868234081,
    
94. 876850801,882796321,885336481,888700681,891706861,897880321,902645857,914801665,
    
95. 918661501,928482241,930745621,931694401,934784929,935794081,939947009,940123801,
    
96. 941056273,945959365,947993761,954732853,955134181,957044881,958735681,958762729,
    
97. 958970545,962442001,962500561,963163201,963168193,968553181,968915521,975303121,
    
98. 977737321,977892241,981567505,981789337,985052881,986088961,990893569,993420289,
    
99. 993905641,1001152801,1018928485,1027334881,1030401901,1031750401,1035608041,
    
100. 1038165961,1055384929,1070659201,1072570801,1074363265,1079556193,1090842145,
     
101. 1093916341,1100674561,1103145121,1125038377,1131222841,1132988545,1134044821,
     
102. 1136739745,1138049137,1140441121,1150270849,1152793621,1162202581,1163659861,
     
103. 1177195201,1177800481,1180398961,1183104001,1189238401,1190790721,1193229577,
     
104. 1194866101,1198650961,1200456577,1200778753,1206057601,1207252621,1210178305,
     
105. 1213619761,1214703721,1216631521,1223475841,1227220801,1227280681,1232469001,
     
106. 1251295501,1251992281,1254318481,1256855041,1257102001,1260332137,1264145401,
     
107. 1268604001,1269295201,1271325841,1295577361,1299963601,1309440001,1312114945,
     
108. 1312332001,1316958721,1317828601,1318126321,1321983937,1330655041,1332521065,
     
109. 1337805505,1348964401,1349671681,1362132541,1376844481,1378483393,1382114881,
     
110. 1384157161,1394746081,1394942473,1404111241,1404228421,1407548341,1410833281,
     
111. 1419339691,1420379065,1422477001,1423668961,1428966001,1439328001,1439492041,
     
112. 1441316269,1442761201,1445084173,1452767521,1481619601,1483199641,1490078305,
     
113. 1490522545,1504651681,1505432881,1507746241,1515785041,1520467201,1521835381,
     
114. 1528936501,1529544961,1534274841,1540454761,1545387481,1571503801,1573132561,
     
115. 1574601601,1576826161,1583582113,1588247851,1592668441,1597009393,1597821121,
     
116. 1615335085,1618206745,1626167341,1632785701,1641323905,1646426881,1648076041,
     
117. 1657700353,1659935761,1672719217,1674309385,1676203201,1676641681,1678569121,
     
118. 1685266561,1688214529,1688639041,1689411601,1690230241,1696572001,1698623641,
     
119. 1699279441,1701016801,1705470481,1708549501,1726372441,1746692641,1750412161,
     
120. 1752710401,1760460481,1769031901,1772267281,1773486001,1776450565,1778382541,
     
121. 1784291041,1784975941,1785507361,1795216501,1801558201,1803278401,1817067169,
     
122. 1825568641,1828377001,1831048561,1833328621,1836304561,1841034961,1845871105,
     
123. 1846817281,1848112761,1848681121,1849811041,1854001513,1855100017,1858395529,
     
124. 1879480513,1887933601,1894344001,1896961801,1899525601,1913016001,1916729101,
     
125. 1918052065,1919767681,1932608161,1942608529,1943951041,1949646601,1950276565,
     
126. 1954174465,1955324449,1958102641,1962804565,1976295241,1984089601,1988071801,
     
127. 1992841201,1999004365,2000436751,2023528501,2029554241,2048443501,2048751901,
     
128. 2049293401,2064236401,2064373921,2067887557,2073560401,2080544005,2097317377,
     
129. 2101170097,2105594401,2107535221,2111416021,2111488561,2117725921,2126689501,
     
130. 2140538401,2140699681,2159003281,2170282969,2176838049,2178944461,2199700321,
     
131. 2199931651,2201169601,2212935985,2215407601,2216430721,2217951073,2223876601,
     
132. 2224519921,2230305949,2232385345,2239622113,2244932281,2246916001,2258118721,
     
133. 2265650401,2272748401,2278677961,2295209281,2301745249,2302419601,2309027281,
     
134. 2313774001,2320224481,2320690177,2322648901,2323147201,2332627249,2335640077,
     
135. 2339165521,2342644921,2353639681,2359686241,2361232477,2367379201,2377166401,
     
136. 2391137281,2396357041,2407376665,2414829781,2430556381,2436691321,2438403661,
     
137. 2443829641,2444950561,2456536681,2457411265,2467813621,2470348441,2470894273,
     
138. 2479305985,2480343553,2489462641,2494465921,2494621585,2497638781,2509860961,
     
139. 2510085721,2519819281,2523947041,2527812001,2529410281,2539024741,2544590161,
     
140. 2547621973,2560104001,2560600351,2561945401,2564889601,2573686441,2574243721,
     
141. 2575260241,2586927553,2588653081,2597928961,2598933481,2601144001,2602378721,
     
142. 2605557781,2607162961,2607237361,2616662881,2617181281,2630374741,2642025673,
     
143. 2657502001,2665141921,2677147201,2685422593,2690867401,2693939401,2702470861,
     
144. 2709611521,2723859001,2733494401,2735309521,2766172501,2766901501,2770560241,
     
145. 2776874941,2787998641,2797002901,2801124001,2806205689,2811315361,2815304401,
     
146. 2832480001,2833846561,2842912381,2858298301,2867755969,2942952481,2943556201,
     
147. 2965085641,2998467901,3001561441,3007991701,3022354401,3024774901,3025708561,
     
148. 3030758401,3034203361,3035837161,3044238121,3044970001,3068534701,3069196417,
     
149. 3072080089,3072094201,3077802001,3078386641,3086434561,3088134721,3090578401,
     
150. 3102234751,3104207821,3105567361,3112974481,3119101921,3138302401,3159422785,
     
151. 3159939601,3164207761,3180288385,3180632833,3188744065,3190894201,3193414093,
     
152. 3203895601,3215031751,3222053185,3232450585,3240392401,3245477761,3246238801,
     
153. 3248891101,3249390145,3263564305,3264820001,3270933121,3277595665,3281736601,
     
154. 3284630713,3307322305,3313196881,3313744561,3314111761,3319323601,3328437481,
     
155. 3345878017,3347570941,3348463105,3353809537,3378014641,3380740301,3411338491,
     
156. 3413656441,3429457921,3438721441,3441837421,3480174001,3504570301,3508507801,
     
157. 3521441665,3534510001,3555636481,3574014445,3575798785,3576804001,3600918181,
     
158. 3618244081,3630291841,3637405045,3682471321,3697952401,3711456001,3712280041,
     
159. 3713287801,3715938721,3722793481,3727589761,3754483201,3767865601,3776698801,
     
160. 3787491457,3799111681,3800513761,3801823441,3805181281,3832646221,3834444901,
     
161. 3835537861,3858853681,3863326897,3880251649,3891338101,3892568065,3901730401,
     
162. 3907357441,3922752121,3928256641,3951813601,3981047941,3998554561,4015029061,
     
163. 4030864201,4034969401,4059151489,4065133501,4077957961,4115677501,4127050621,
     
164. 4134273793,4138747921,4146685921,4160472121,4162880401,4167038161,4169092201,
     
165. 4169867689,4189909501,4199202001,4199529601,4199932801,4202009461,4210922233,
     
166. 4212413569,4215885697,4216799521,4277982241        // 1119: 4295605861, 溢出
     
167. };
     
168. 
     
169. // hash表的改进方法, 只需要保存下标
     
170. // 经测试, 最多有6次冲突, b[i][0]保存数目
     
171. // 注意buckets的各维大小是人工测试得到!!
     
172. 
     
173. static short buckets[1021][6+1];
     
174. 
     
175. // 自动初始化hash表
     
176. 
     
177. static struct init
     
178. {
     
179.         init(void)
     
180.         {
     
181.                 for(int i = 0; i < NELEMS(CarmichaeNumbers); ++i)
     
182.                 {
     
183.                         // buckets[v]类似一个栈
     
184.                         // buckets[v][0]记录栈顶位置
     
185.                         // 最多有6次冲突, 不会溢出
     
186. 
     
187.                         int key = CarmichaeNumbers[i]%NELEMS(buckets);
     
188.                         assert(buckets[key][0] < NELEMS(buckets[0])-1);
     
189. 
     
190.                         buckets[key][++buckets[key][0]] = i;
     
191.                 }
     
192.         }
     
193. } _init;
     
194. 
     
195. // 判断是否是卡尔麦克数
     
196. 
     
197. static bool isCarmichaelNum(unsigned n)
     
198. {
     
199.         // 改进: 采用hash表查找
     
200. 
     
201.         int i, key = n%NELEMS(buckets);
     
202. 
     
203.         for(i = 1; i <= buckets[key][0]; ++i)
     
204.                 if(n == CarmichaeNumbers[buckets[key][i]]) return true;
     
205. 
     
206.         return false;
     
207. }
     
208. 
     
209. // 计算(u*v)%m
     
210. 
     
211. unsigned mul_mod(unsigned u, unsigned v, unsigned z)
     
212. {
     
213.         // 如果u*v没有溢出, 则直接计算
     
214. 
     
215.         if((u*v)/u == v) return (u*v)%z;
     
216. 
     
217.         // 进行长乘法(结果为64位)
     
218. 
     
219.         unsigned u0, v0, w0;
     
220.         unsigned u1, v1, w1, w2, t;
     
221. 
     
222.         u0 = u & 0xFFFF;  u1 = u >> 16;
     
223.         v0 = v & 0xFFFF;  v1 = v >> 16;
     
224.         w0 = u0*v0;
     
225.         t  = u1*v0 + (w0 >> 16);
     
226.         w1 = t & 0xFFFF;
     
227.         w2 = t >> 16;
     
228.         w1 = u0*v1 + w1;
     
229. 
     
230.         // x为高32位, y为低32位
     
231. 
     
232.         unsigned x = u1*v1 + w2 + (w1 >> 16);
     
233.         unsigned y = u*v;
     
234. 
     
235.         // 进行长除法(被除数64位, 除数32位)
     
236. 
     
237.         for (int i = 1; i <= 32; i++)
     
238.         {
     
239.                 t = (int)x >> 31;           // All 1's if x(31) = 1.
     
240. 
     
241.                 x = (x << 1) | (y >> 31);   // Shift x || y left
     
242.                 y <<= 1;                    // one bit.
     
243.                
     
244.                 if((x|t) >= z) { x -= z; y++; }
     
245.         }
     
246. 
     
247.         return x;        // y为商, x为余数
     
248. }
     
249. 
     
250. // 判断是否是费马(伪)素数
     
251. 
     
252. static bool isFermatNum(unsigned p)
     
253. {
     
254.         unsigned m = p--;
     
255.         unsigned r = 2%m;
     
256.         unsigned k = 1;
     
257.    
     
258.         // 蒙格马利快速幂模算法
     
259. 
     
260.         while(p > 1)
     
261.         {
     
262.                 if(p&1) k = mul_mod(k, r, m);
     
263.                 r = mul_mod(r, r, m); p >>= 1;
     
264.         }
     
265. 
     
266.         return mul_mod(k, r, m) == 1;
     
267. }
     
268. 
     
269. // 根据费马小定理测试
     
270. 
     
271. bool isPrime(unsigned n)
     
272. {
     
273.         // 处理比较特殊的情况
     
274. 
     
275.         if(n < 2) return false;
     
276.         if(n == 2) return true;
     
277.         if(!(n&1)) return false;
     
278. 
     
279.         // 是费马(伪)素数但不是卡尔麦克数则必定为素数
     
280. 
     
281.         return isFermatNum(n) && !isCarmichaelNum(n);
     
282. }
     

复制代码

这个方法其实还有缺陷，主要因为之前概念理解有错误。
卡尔麦克数是针对所有Zn*费马测试为真的情况！
因此，基于2（或者其他值）的伪素数，除了卡尔麦克数之外还有别的。
比如：341。这个之前都没有考虑到。

不过，即使只针对2的情况，伪素数还是比较少的，
可以在我们期望的一个范围内。
那样，依然可以保存基于2的所有伪素数表。

当然，生成2^32之内的2的伪素数可能要花些时间
（我采用的是比较普通的方法）

下次我将给出生成2的伪素数和一些自动测试的代码

chai2010

开始的方法代码也可以这样

3. // 计算数组的元素个数
5. #define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))
7. // 构造素数序列primes[]
9. void makePrimes(int primes[], int num)
10. {
11.         primes[0] = 2;
12.    
13.         for(int p = 3, cnt = 1; cnt < num; p += 2)
14.         {
15.                 for(int k = 1; ; ++k)
16.                 {
17.                         // 可以在div时得到余数
18.                         // assert(p <= INT_MAX);
20.                         div_t dt = div(p, primes[k]);
22.                         if(dt.rem == 0) break;
23.                         if(dt.quot <= primes[k]) { primes[cnt++] = p; break; }
24.                 }
25.         }
26. }
28. // 其他算法判断素数, 用于测试
30. bool isPrime(unsigned p)
31. {
32.         // 静态初始化素数序列
34.         static int px[1024*8];
35.         static struct Init {
36.                 Init(int *p, int n){makePrimes(p, n);}
37.         } _init(px, NELEMS(px));
39.         // 处理特殊情况
41.         if(p < 2) return false;
42.         if(p == 2) return true;
43.         if(!(p&1)) return false;
45.         // 测试px[i]是否为其因子
47.         for(int k = 1; px[k]*px[k] <= p; ++k)
48.         {
49.                 if(p%px[k] == 0) return false;
50.         }
51.         return true;
52. }

复制代码

chai2010

测试程序

1. 
   
2. #include <assert.h>
   
3. #include <stdlib.h>
   
4. #include <time.h>
   
5. 
   
6. bool Test_isPrime(unsigned p)
   
7. {
   
8.         // 其他方法判断
   
9. }
   
10. 
    
11. int main(int max)
    
12. {
    
13.         srand(time(NULL));
    
14. 
    
15.         for(unsigned i = 0; i < max; ++i)
    
16.         {
    
17.                 unsigned p = (rand()<<16) | rand();
    
18.                 assert(isPrime(p) == Test_isPrime(p));
    
19.         }
    
20. }
    

复制代码

chai2010

生成2^32以内基于2的费马伪素数

1. 
   
2. // 缺少的函数参考前面的帖子
   
3. // 测试(1, 2^32)内所有基于2的费马伪素数
   
4. // 偶数的情况没有测试
   
5. 
   
6. int main()
   
7. {
   
8.         unsigned i, cnt = 0;
   
9. 
   
10.         for(i = 1; i <= UINT_MAX; i += 2)
    
11.         {
    
12.                 if(isPrime(i) && (!Test_isPrime(i)))
    
13.                 {
    
14.                         printf("%d: %u\n", ++cnt, i);
    
15.                 }
    
16.         }
    
17. 
    
18.         return 0;
    
19. }

复制代码

可以将输出重定向到文件中: ./a.out > data.txt

由于耗时比较长，最好是在机器空闲的时候跑:-)
我打算晚上让机器跑一夜，有结果在发上来...

chai2010

这里有基于2的费马伪素数，

http://www.research.att.com/~njas/sequences/b001567.txt

chai2010

修改后的代码

1. 
   
2. 
   
3. // 函数: bool isPrime(unsigned n);
   
4. // 描述: 判断n(0, 2^32-1)是否为素数
   
5. // 算法: 费马小定理剔除伪素数
   
6. // 作者: 柴树杉
   
7. // 时间: 2006-11-06
   
8. // 主页: [url]chaishushan.googlepages.com[/url]
   
9. 
   
10. #include "FermatPrime.h"
    
11. 
    
12. #include <assert.h>
    
13. #include <stdlib.h>
    
14. 
    
15. // 计算数组的元素个数
    
16. 
    
17. #define NELEMS(x) ((sizeof(x)) / (sizeof((x)[0])))
    
18. 
    
19. // 2^32之内基于2的费马伪素数
    
20. 
    
21. static const unsigned FermatPseudoprimes[] =
    
22. {
    
23. 341,561,645,1105,1387,1729,1905,2047,2465,2701,2821,3277,4033,4369,4371,
    
24. 4681,5461,6601,7957,8321,8481,8911,10261,10585,11305,12801,13741,13747,
    
25. 13981,14491,15709,15841,16705,18705,18721,19951,23001,23377,25761,29341,
    
26. 30121,30889,31417,31609,31621,33153,34945,35333,39865,41041,41665,42799,
    
27. 46657,49141,49981,52633,55245,57421,60701,60787,62745,63973,65077,65281,
    
28. 68101,72885,74665,75361,80581,83333,83665,85489,87249,88357,88561,90751,
    
29. 91001,93961,101101,104653,107185,113201,115921,121465,123251,126217,129889,
    
30. 129921,130561,137149,149281,150851,154101,157641,158369,162193,162401,164737,
    
31. 
    
32. // ... 缺少的伪素数请访问:
    
33. // [url]http://www.research.att.com/~njas/sequences/b001567.txt[/url]
    
34. 
    
35. 4284050473,4285148981,4286383201,4286813749,4288664869,4289470021,4289641621,
    
36. 4289884201,4289906089,4293088801,4293329041,4294868509,4294901761
    
37. };        // num = 10403
    
38. 
    
39. // 构造hash表
    
40. 
    
41. static short buckets[10427];
    
42. static short link[NELEMS(FermatPseudoprimes)];
    
43. 
    
44. // 建立hash表, 只运行一次
    
45. // hash的最大冲突为7次
    
46. 
    
47. static void initBuckets(void)
    
48. {
    
49.         int max = 0;
    
50.         for(int i = 0; i < NELEMS(FermatPseudoprimes); ++i)
    
51.         {
    
52.                 int h = FermatPseudoprimes[i]%NELEMS(buckets);
    
53. 
    
54.                 { link[i] = buckets[h]; buckets[h] = i + 1; }
    
55.         }
    
56. }
    
57. 
    
58. // 判断是否是费马伪素数
    
59. 
    
60. static bool isFermatPseudoprime(unsigned p)
    
61. {
    
62.         // 静态调用initBuckets函数
    
63. 
    
64.         static struct Init {
    
65.                 Init(){initBuckets();} } _init;
    
66. 
    
67.         // 在hash表中查找
    
68. 
    
69.         int h = p%NELEMS(buckets);
    
70.         int next = buckets[h]-1;
    
71. 
    
72.         for( ; next >= 0; next = link[next]-1)
    
73.         {
    
74.                 if(p == FermatPseudoprimes[next]) return true;
    
75.         }
    
76.         return false;
    
77. }
    
78. 
    
79. // 计算(u*v)%m
    
80. 
    
81. unsigned mul_mod(unsigned u, unsigned v, unsigned z)
    
82. {
    
83.         // 如果u*v没有溢出, 则直接计算
    
84. 
    
85.         if((u*v)/u == v) return (u*v)%z;
    
86. 
    
87.         // 进行长乘法(结果为64位)
    
88. 
    
89.         unsigned u0, v0, w0;
    
90.         unsigned u1, v1, w1, w2, t;
    
91. 
    
92.         u0 = u & 0xFFFF;  u1 = u >> 16;
    
93.         v0 = v & 0xFFFF;  v1 = v >> 16;
    
94.         w0 = u0*v0;
    
95.         t  = u1*v0 + (w0 >> 16);
    
96.         w1 = t & 0xFFFF;
    
97.         w2 = t >> 16;
    
98.         w1 = u0*v1 + w1;
    
99. 
    
100.         // x为高32位, y为低32位
     
101. 
     
102.         unsigned x = u1*v1 + w2 + (w1 >> 16);
     
103.         unsigned y = u*v;
     
104. 
     
105.         // 进行长除法(除数为64位)
     
106. 
     
107.         for (int i = 1; i <= 32; i++)
     
108.         {
     
109.                 t = (int)x >> 31;           // All 1's if x(31) = 1.
     
110. 
     
111.                 x = (x << 1) | (y >> 31);   // Shift x || y left
     
112.                 y <<= 1;                    // one bit.
     
113.                
     
114.                 if((x|t) >= z) { x -= z; y++; }
     
115.         }
     
116. 
     
117.         return x;        // y为商, x为余数
     
118. }
     
119. 
     
120. // 判断是否是费马(伪)素数
     
121. 
     
122. static bool isFermatNum(unsigned p)
     
123. {
     
124.         unsigned m = p--;
     
125.         unsigned r = 2%m;
     
126.         unsigned k = 1;
     
127.    
     
128.         // 蒙格马利快速幂模算法
     
129. 
     
130.         while(p > 1)
     
131.         {
     
132.                 if(p&1) k = mul_mod(k, r, m);
     
133.                 r = mul_mod(r, r, m); p >>= 1;
     
134.         }
     
135. 
     
136.         return mul_mod(k, r, m) == 1;
     
137. }
     
138. 
     
139. // 根据费马小定理测试
     
140. 
     
141. bool isPrime(unsigned n)
     
142. {
     
143.         // 处理比较特殊的情况
     
144. 
     
145.         if(n < 2) return false;
     
146.         if(n == 2) return true;
     
147.         if(!(n&1)) return false;
     
148. 
     
149.         // 从费马(伪)素数剔除伪素数情况
     
150. 
     
151.         return isFermatNum(n) && !isFermatPseudoprime(n);
     
152. }
     
153. 
     
154. 
     

复制代码

chai2010

在第1帖 9 中， 作为副产品曾提供以下2个函数：

extern int  Prime_max(void);        // 素数序列的大小
extern int  Prime_get (int i);        // 返回第i个素数, 0 <= i < Prime_max

Prime_max返回缓冲的素数序列大小，
Prime_get则可以从缓冲的序列中获取一定的素数。

其实Prime_get完全可以返回2^32范围内任意素数！

另外，还可以提供Prime_num(unsigned max);用于返回
1-max之间素数的个数（也用于实现Prime_get函数）。

具体的细节和代码下次慢慢补充...