估计下找到10^100以上的一个四生素数的工作量

gxqcn

原帖由 无心人 于 2008-7-16 10:28 发表


GxQ尝试过多大的数字的素性测试
具体时间是多少？

没具体测试 HugeCalc 素性判定的能力，
应该说是足够大的，可以容纳当前已知的最大素数。

其算法还可进一步优化，留待下一版吧。。。

无心人

我就想知道进行一个10000位数字的一次测试需要多少秒？
是一次测试，即单独一个素数的测试

gxqcn

HugeCalc 内部针对不同的数据特点有不同的加速算法，
比如判定 $10^10000 +- \delta \  ( \delta\ "is a small integer")$ 将比普通的 10000 位整数快很多。

无心人

那就按判定普通数字测试下

无心人

我评估了下素性测试的时间
因为是linux系统，用的GMP

Input:67
Begin Test 2^67 - 1.
Composite Number.
Used Time:0.000097
Input:31  
Begin Test 2^31 - 1.
Probably prime Number.
Used Time:0.000086
Input:127
Begin Test 2^127 - 1.
Probably prime Number.
Used Time:0.000119
Input:61
Begin Test 2^61 - 1.
Probably prime Number.
Used Time:0.000070
Input:257
Begin Test 2^257 - 1.
Composite Number.
Used Time:0.000249
Input:4423
Begin Test 2^4423 - 1.
Probably prime Number.
Used Time:1.412644
Input:9941
Begin Test 2^9941 - 1.
Probably prime Number.
Used Time:11.602822
Input:9689
Begin Test 2^9689 - 1.
Probably prime Number.
Used Time:10.859898
Input:0

无心人

可以看到平均每10秒能测试出一个3000位十进制整数的素性

现在问题是3000位的四生素数
如果用
a * b# + C, a = 1..999999999999, b = 5000..9999, C是所有10^7以下的四生素数的集合
请问，如果以10秒做一次素性筛选
大概多长时间能找到一组通过素性筛选的四生素数可能组合？

gxqcn

GMP 采用的是5次随机 Miller-Rabin 测试，强度并不够，
HugeCalc 采用了更高强度的素性检测。

注：65# 的测试数据都比较特殊，HugeCalc 可以自动切换高速算法检测。

无心人

我不要强度，我要速度
成万的测试，速度第一

另外，不要考虑特殊形式
因为就是为了测试普通数字

无心人

带状态保存的搜索大四生素数程序
最大到(2^31-1) * 9999# + 10^8

1. #include <stdlib.h>
   
2. #include <stdio.h>
   
3. #include <string.h>
   
4. #include <math.h>
   
5. #include <assert.h>
   
6. #include <sys/time.h>
   
7. #include "gmp.h"
   
8. 
   
9. unsigned long primes[1280];        //小于10000素数，实际上有1229个
   
10. unsigned long quadprimes[4800];        //小于10^8四生素数，实际上有4768组
    
11. unsigned long sieve[30030];
    
12. unsigned long template[30030];
    
13. unsigned long A, B, C;
    
14. 
    
15. void
    
16. InitPrimes (void)
    
17. {
    
18.   primes[0] = 2;
    
19.   primes[1] = 3;
    
20.   primes[2] = 5;
    
21.   primes[3] = 7;
    
22.   primes[4] = 11;
    
23.   primes[5] = 13;
    
24.   primes[6] = 17;
    
25.   primes[7] = 19;
    
26.   primes[8] = 23;
    
27.   primes[9] = 29;
    
28.   primes[10] = 31;
    
29.   primes[11] = 37;
    
30.   primes[12] = 41;
    
31.   primes[13] = 43;
    
32.   primes[14] = 47;
    
33.   primes[15] = 53;
    
34.   primes[16] = 59;
    
35.   primes[17] = 61;
    
36.   primes[18] = 67;
    
37.   primes[19] = 71;
    
38.   primes[20] = 73;
    
39.   primes[21] = 79;
    
40.   primes[22] = 83;
    
41.   primes[23] = 89;
    
42.   primes[24] = 97;
    
43. }
    
44. 
    
45. void
    
46. InitSieve (void)
    
47. {
    
48.   unsigned long i;
    
49.   for (i = 1; i <= 9999; i += 2)
    
50.     sieve[i] = 1;
    
51. }
    
52. 
    
53. void
    
54. Sieve (void)
    
55. {
    
56.   unsigned long i, j, start;
    
57.   for (i = 0; i <= 24; i++)
    
58.     {
    
59.       start = primes[i] * primes[i];
    
60.       for (j = start; j <= 9999; j += 2 * primes[i])
    
61.         sieve[j] = 0;
    
62.     }
    
63. }
    
64. 
    
65. void
    
66. StoreSieveResult (void)
    
67. {
    
68.   unsigned long i, k;
    
69.   k = 25;
    
70.   for (i = 101; i <= 9999; i += 2)
    
71.     if (sieve[i])
    
72.       primes[k++] = i;
    
73.   for (; k < 1280; k++)
    
74.     primes[k] = 0;
    
75.   primes[1229] = 10007;                //以这个素数结束，否则下面的四生素数筛在接近10^8会超越1229素数的范围
    
76. }
    
77. 
    
78. void
    
79. InitQuadPrimes (void)
    
80. {
    
81.   unsigned long i, k = 0;
    
82.   for (i = 0; i < 1229; i++)
    
83.     if ((primes[i + 1] - primes[i] == 2) && (primes[i + 2] - primes[i] == 6)
    
84.         && (primes[i + 3] - primes[i] == 8))
    
85.       quadprimes[k++] = primes[i];        //共12组
    
86. }
    
87. 
    
88. void
    
89. InitTemplate (void)
    
90. {
    
91.   unsigned long i;
    
92.   for (i = 1; i < 30030; i += 2)
    
93.     template[i] = 1;
    
94.   for (i = 0; i < 30030; i += 2)
    
95.     template[i] = 0;
    
96.   for (i = 3; i < 30030; i += 6)
    
97.     template[i] = 0;
    
98.   for (i = 5; i < 30030; i += 10)
    
99.     template[i] = 0;
    
100.   for (i = 7; i < 30030; i += 14)
     
101.     template[i] = 0;
     
102.   for (i = 11; i < 30030; i += 22)
     
103.     template[i] = 0;
     
104.   for (i = 13; i < 30030; i += 26)
     
105.     template[i] = 0;
     
106. }
     
107. 
     
108. void
     
109. SearchQuadSieve (void)
     
110. {
     
111.   unsigned long start, end;
     
112.   unsigned long i, j, k, m, s_start;
     
113. 
     
114.   InitTemplate ();
     
115. 
     
116. //首次筛过程
     
117.   start = 0;
     
118.   end = 30030;
     
119.   memcpy (sieve, template, 30030 * sizeof (unsigned long));
     
120. 
     
121.   for (j = 6; primes[j] * primes[j] <= end; j++)
     
122.     {
     
123.       s_start = primes[j] * primes[j];
     
124.       for (k = s_start; k < 30030; k += 2 * primes[j])
     
125.         sieve[k] = 0;
     
126.     }
     
127. 
     
128.   m = 12;
     
129.   for (i = 10001; i < 30030; i += 10)
     
130.     if ((sieve[i]) && (sieve[i + 2]) && (sieve[i + 6]) && (sieve[i + 8]))
     
131.       quadprimes[m++] = i;
     
132. 
     
133. //完整筛过程，经测试剩余的数字是99999900-99999999，不存在四生素数
     
134.   for (i = 1; i < 100000000 / 30030; i++)
     
135.     {
     
136.       start = i * 30030;
     
137.       end = start + 30030;
     
138.       memcpy (sieve, template, 30030 * sizeof (unsigned long));
     
139. 
     
140.       for (j = 6; primes[j] * primes[j] <= start; j++)
     
141.         {
     
142.           s_start = (start / primes[j]) * primes[j] + primes[j];
     
143.           if (s_start % 2 == 0)
     
144.             s_start += primes[j];
     
145.           s_start -= start;
     
146.           for (k = s_start; k < 30030; k += 2 * primes[j])
     
147.             sieve[k] = 0;
     
148.         }
     
149. 
     
150.       for (; primes[j] * primes[j] < end; j++)
     
151.         {
     
152.           s_start = primes[j] * primes[j] - start;
     
153.           for (k = s_start; k < 30030; k += 2 * primes[j])
     
154.             sieve[k] = 0;
     
155.         }
     
156. 
     
157.       s_start = (start / 30) * 30 + 11;
     
158.       if (s_start < start)
     
159.         s_start += 30;
     
160.       s_start -= start;
     
161. 
     
162.       for (k = s_start; k < 30030; k += 30)
     
163.         if ((sieve[k]) && (sieve[k + 2]) && (sieve[k + 6]) && (sieve[k + 8]))
     
164.           quadprimes[m++] = start + k;
     
165. 
     
166.     }
     
167. }
     
168. 
     
169. void
     
170. Init (void)
     
171. {
     
172.   struct timeval start, end;
     
173.   double timeuse;
     
174. 
     
175.   gettimeofday (&start, NULL);
     
176.   InitPrimes ();
     
177.   InitSieve ();
     
178.   Sieve ();
     
179.   StoreSieveResult ();
     
180. 
     
181.   InitQuadPrimes ();
     
182.   SearchQuadSieve ();
     
183.   gettimeofday (&end, NULL);
     
184. 
     
185.   timeuse =
     
186.     1000000 * (end.tv_sec - start.tv_sec) + end.tv_usec - start.tv_usec;
     
187.   timeuse /= 1000000;
     
188. 
     
189.   printf ("Time use:%8.6f:  %\n", timeuse);
     
190. 
     
191. }
     
192. 
     
193. void
     
194. product (mpz_t r, unsigned long p[], unsigned long max)
     
195. {
     
196.   unsigned long i;
     
197.   mpz_set_ui (r, 1);
     
198. 
     
199.   if (max < 2)
     
200.     return;
     
201. 
     
202.   for (i = 0; primes[i] <= max; i++)
     
203.     mpz_mul_ui (r, r, p[i]);
     
204. }
     
205. 
     
206. //状态文件分三行，分别存储A, B, C的数值
     
207. void
     
208. LoadStatus (void)
     
209. {
     
210.   FILE *fstatus;
     
211.   fstatus = fopen ("quadprimes.status", "r+t");
     
212.   if (fstatus == NULL)
     
213.     {
     
214.       printf ("读取状态文件失败。\n");
     
215.       exit (1);
     
216.     }
     
217.   fscanf (fstatus, "%u%u%u", &A, &B, &C);
     
218.   fclose (fstatus);
     
219. }
     
220. 
     
221. void
     
222. SaveStatus (void)
     
223. {
     
224.   FILE *fstatus;
     
225.   fstatus = fopen ("quadprimes.status", "w+t");
     
226.   if (fstatus == NULL)
     
227.     {
     
228.       printf ("存储状态文件失败。\n");
     
229.       exit (1);
     
230.     }
     
231.   fprintf (fstatus, "%u\n%u\n%u\n", A, B, C);
     
232.   fclose (fstatus);
     
233. }
     
234. 
     
235. void
     
236. SaveLog (void)
     
237. {
     
238.   FILE *flog;
     
239.   flog = fopen ("quadprimes.log", "a+t");
     
240.   if (flog == NULL)
     
241.     {
     
242.       printf ("存储日志文件失败。\n");
     
243.       exit (1);
     
244.     }
     
245.   fprintf (flog, "%u    %u    %u\n", A, B, C);
     
246.   printf ("%u    %u    %u\n", A, B, C);
     
247.   fclose (flog);
     
248. 
     
249. }
     
250. 
     
251. unsigned long SearchQuadPrime(unsigned long p)
     
252. {
     
253.   unsigned long k;
     
254.   for (k = 0; k < 4768; k ++)
     
255.      if (quadprimes[k] >= p)
     
256.           return (k);
     
257.     return (0);
     
258. }
     
259. 
     
260. void
     
261. Test (void)
     
262. {
     
263.   unsigned long m, k;
     
264.   mpz_t b0, b, p1, p2, p6, p8;
     
265.   LoadStatus ();
     
266.   mpz_init (b0);
     
267.   mpz_init(b);
     
268.   mpz_init (p1);
     
269.   mpz_init (p2);
     
270.   mpz_init (p6);
     
271.   mpz_init (p8);
     
272. 
     
273.   product (b0, primes, B);
     
274.   k = SearchQuadPrime(C);
     
275.   C = quadprimes[k];
     
276.   mpz_mul_ui (b, b0, A);
     
277.   mpz_add_ui (p1, b, C);
     
278.   m = 0;
     
279.   for (;;)
     
280.     {
     
281.       if (mpz_probab_prime_p (p1, 1))
     
282.         {
     
283.           mpz_add_ui (p2, p1, 2);
     
284.           if (mpz_probab_prime_p (p2, 1))
     
285.             {
     
286.               mpz_add_ui (p6, p1, 6);
     
287.               if (mpz_probab_prime_p (p6, 1))
     
288.                 {
     
289.                   mpz_add_ui (p8, p1, 8);
     
290.                   if (mpz_probab_prime_p (p8, 1))
     
291.                  SaveLog ();
     
292.                 }
     
293.             }
     
294.          }
     
295. 
     
296.     //
     
297.         m ++;
     
298.         if (m > 100)
     
299.           {
     
300.           SaveStatus ();
     
301.           m = 0;
     
302.           }
     
303.         k ++;
     
304.         if (k >= 4768)
     
305.           {
     
306.           k = SearchQuadPrime(B);
     
307.           C = quadprimes[k];
     
308.           A ++;
     
309.           if (A == 0)
     
310.              {
     
311.             SaveLog();
     
312.             exit(2);
     
313.              }
     
314.           else
     
315.              {
     
316.             mpz_add(b, b, b0);
     
317.             mpz_add_ui(p1, b, C);
     
318.              }
     
319.           }
     
320.         else
     
321.           {
     
322.           C =  quadprimes[k];
     
323.           mpz_add_ui(p1, b, C);
     
324.           }
     
325.     }
     
326.   mpz_clear (b);
     
327.   mpz_clear(b0);
     
328.   mpz_clear (p1);
     
329.   mpz_clear (p2);
     
330.   mpz_clear (p6);
     
331.   mpz_clear (p8);
     
332. }
     
333. 
     
334. int main (void)
     
335. {
     
336.   unsigned long a, c;
     
337. 
     
338.   Init ();
     
339. 
     
340.   Test ();
     
341. //  for (a = 0; a < 4800; a++)
     
342. //    if (quadprimes[a])
     
343. //      printf ("%6u ", quadprimes[a]);
     
344. //  printf ("\n");
     
345.   return 0;
     
346. }

复制代码

无心人

上面这个程序
我忘记是否是调通的了

总想拿这个来算下