孪生素数的计算

无心人

谢谢

我在写一个东西
准备挂服务器上
算2000－3000位的四生素数
呵呵

无心人

1. #include <stdlib.h>
   
2. #include <stdio.h>
   
3. #include <string.h>
   
4. #include <math.h>
   
5. #include <assert.h>
   
6. #include "gmp.h"
   
7. 
   
8. unsigned long primes[1280];        //小于10000素数，实际上有1229个
   
9. unsigned long quadprimes[4800];        //小于10^8四生素数，实际上有4768组
   
10. unsigned long sieve[30030];
    
11. unsigned long template[30030];
    
12. 
    
13. void
    
14. InitPrimes (void)
    
15. {
    
16.   primes[0] = 2;
    
17.   primes[1] = 3;
    
18.   primes[2] = 5;
    
19.   primes[3] = 7;
    
20.   primes[4] = 11;
    
21.   primes[5] = 13;
    
22.   primes[6] = 17;
    
23.   primes[7] = 19;
    
24.   primes[8] = 23;
    
25.   primes[9] = 29;
    
26.   primes[10] = 31;
    
27.   primes[11] = 37;
    
28.   primes[12] = 41;
    
29.   primes[13] = 43;
    
30.   primes[14] = 47;
    
31.   primes[15] = 53;
    
32.   primes[16] = 59;
    
33.   primes[17] = 61;
    
34.   primes[18] = 67;
    
35.   primes[19] = 71;
    
36.   primes[20] = 73;
    
37.   primes[21] = 79;
    
38.   primes[22] = 83;
    
39.   primes[23] = 89;
    
40.   primes[24] = 97;
    
41. }
    
42. 
    
43. void
    
44. InitSieve (void)
    
45. {
    
46.   unsigned long i;
    
47.   for (i = 1; i <= 9999; i += 2)
    
48.     sieve[i] = 1;
    
49. }
    
50. 
    
51. void
    
52. Sieve (void)
    
53. {
    
54.   unsigned long i, j, start;
    
55.   for (i = 0; i <= 24; i++)
    
56.     {
    
57.       start = primes[i] * primes[i];
    
58.       for (j = start; j <= 9999; j += 2 * primes[i])
    
59.         sieve[j] = 0;
    
60.     }
    
61. }
    
62. 
    
63. void
    
64. StoreSieveResult (void)
    
65. {
    
66.   unsigned long i, k;
    
67.   k = 25;
    
68.   for (i = 101; i <= 9999; i += 2)
    
69.     if (sieve[i])
    
70.       primes[k++] = i;
    
71.   for (; k < 1280; k++)
    
72.     primes[k] = 0;
    
73.   primes[1229] = 10007;                //以这个素数结束，否则下面的四生素数筛在接近10^8会超越1229素数的范围
    
74. }
    
75. 
    
76. void
    
77. Init (void)
    
78. {
    
79.   InitPrimes ();
    
80.   InitSieve ();
    
81.   Sieve ();
    
82.   StoreSieveResult ();
    
83. }
    
84. 
    
85. void
    
86. InitQuadPrimes (void)
    
87. {
    
88.   unsigned long i, k = 0;
    
89.   for (i = 0; i < 1229; i++)
    
90.     if ((primes[i + 1] - primes[i] == 2) && (primes[i + 2] - primes[i] == 6)
    
91.         && (primes[i + 3] - primes[i] == 8))
    
92.       quadprimes[k++] = primes[i];        //共12组
    
93. }
    
94. 
    
95. void
    
96. InitTemplate (void)
    
97. {
    
98.   unsigned long i;
    
99.   for (i = 1; i < 30030; i += 2)
    
100.     template[i] = 1;
     
101.   for (i = 0; i < 30030; i += 2)
     
102.     template[i] = 0;
     
103.   for (i = 3; i < 30030; i += 6)
     
104.     template[i] = 0;
     
105.   for (i = 5; i < 30030; i += 10)
     
106.     template[i] = 0;
     
107.   for (i = 7; i < 30030; i += 14)
     
108.     template[i] = 0;
     
109.   for (i = 11; i < 30030; i += 22)
     
110.     template[i] = 0;
     
111.   for (i = 13; i < 30030; i += 26)
     
112.     template[i] = 0;
     
113. }
     
114. 
     
115. void
     
116. SearchQuadSieve (void)
     
117. {
     
118.   unsigned long start, end;
     
119.   unsigned long i, j, k, m, s_start;
     
120. 
     
121.   InitTemplate ();
     
122. 
     
123. //首次筛过程
     
124.   start = 0;
     
125.   end = 30030;
     
126.   memcpy (sieve, template, 30030 * sizeof (unsigned long));
     
127. 
     
128.   for (j = 6; primes[j] * primes[j] <= end; j++)
     
129.     {
     
130.       s_start = primes[j] * primes[j];
     
131.       for (k = s_start; k < 30030; k += 2 * primes[j])
     
132.         sieve[k] = 0;
     
133.     }
     
134. 
     
135.   m = 12;
     
136.   for (i = 10001; i < 30030; i += 10)
     
137.     if ((sieve[i]) && (sieve[i + 2]) && (sieve[i + 6]) && (sieve[i + 8]))
     
138.       quadprimes[m++] = i;
     
139. 
     
140. //完整筛过程，经测试剩余的数字是99999900-99999999，不存在四生素数
     
141.   for (i = 1; i < 100000000 / 30030; i++)
     
142.     {
     
143.       start = i * 30030;
     
144.       end = start + 30030;
     
145.       memcpy (sieve, template, 30030 * sizeof (unsigned long));
     
146. 
     
147.       for (j = 6; primes[j] * primes[j] <= start; j++)
     
148.         {
     
149.           s_start = (start / primes[j]) * primes[j] + primes[j];
     
150.           if (s_start % 2 == 0)
     
151.             s_start += primes[j];
     
152.           s_start -= start;
     
153.           for (k = s_start; k < 30030; k += 2 * primes[j])
     
154.             sieve[k] = 0;
     
155.         }
     
156. 
     
157.       for (; primes[j] * primes[j] < end; j++)
     
158.         {
     
159.           s_start = primes[j] * primes[j] - start;
     
160.           for (k = s_start; k < 30030; k += 2 * primes[j])
     
161.             sieve[k] = 0;
     
162.         }
     
163. 
     
164.       s_start = (start / 30) * 30 + 11;
     
165.       if (s_start < start)
     
166.         s_start += 30;
     
167.       s_start -= start;
     
168. 
     
169.       for (k = s_start; k < 30030; k += 30)
     
170.         if ((sieve[k]) && (sieve[k + 2]) && (sieve[k + 6]) && (sieve[k + 8]))
     
171.           quadprimes[m++] = start + k;
     
172. 
     
173.     }
     
174. }
     
175. 
     
176. void
     
177. product (mpz_t r, unsigned long p[], unsigned long max)
     
178. {
     
179.   unsigned long i;
     
180.   mpz_set_ui (r, 1);
     
181. 
     
182.   if (max < 2)
     
183.     return;
     
184. 
     
185.   for (i = 0; primes[i] <= max; i++)
     
186.     mpz_mul_ui (r, r, p[i]);
     
187. }
     
188. 
     
189. int
     
190. main (void)
     
191. {
     
192.   unsigned long a, c;
     
193.   mpz_t b, p;
     
194. 
     
195.   Init ();
     
196. 
     
197.   InitQuadPrimes ();
     
198.   SearchQuadSieve ();
     
199. 
     
200.   for (a = 0; a < 4800; a++)
     
201.     if (quadprimes[a])
     
202.       printf ("%6u ", quadprimes[a]);
     
203.   printf ("\n");
     
204. 
     
205.   return 0;
     
206. }

复制代码

无心人

上面是俺写大四生素数到初始化部分的程序
初始化了一个10000内素数表，一个100000000内四生素数表
下面是四生素数的结果，不知道对不对

quadprimes.txt.bz2 (16.21 KB, 下载次数: 4)

2008-7-17 10:29 上传

点击文件名下载附件

无心人

不包括输出是570ms
P4 2.0G/1280M/ubuntu 8.04/GCC 4.2.3

tprime

这里有10亿的

tprime

眼馋你的服务器, 多少个节点? CPU配置?

无心人

一个破服务器啊
XEON 2.0 *2/1024M
就是可24小时开机而已

谢谢你的文档

tprime

终于完成计算[0-2^32]任意两数之间的孪生素数(比最快的版本慢1倍, 不同算法)
随机测试100万组数据还没发现错误.目前还不能支持孪生素数列表(输出太慢不打算做了)
欢迎大家下载玩并提出意见, 下面是截屏输出和自己的测试数据.

sieve cache size = 2 * 30 * 7 * 11 * 13 * 17 * 12 / 60 = 199K
PI2(2147483647) = 6810670, init twin prime table time use 1767.14 ms
please input two number: (0<= beg <= end < 2^32):

1 2147483647
case 1 : PI2[1,2147483647] = 6810670, time use 3 ms
1234567890 2147483647
case 2 : PI2[1234567890,2147483647] = 2672471, time use 6 ms
123456789 987654321
case 3 : PI2[123456789,987654321] = 2856051, time use 5 ms

独木星空

8楼的数据很有用

独木星空

N值→→→→4(8)素数个数→孪生素数对数目→2.38128115124→误差→误差率CP
10000000→898→→→→→58755→→→→822→→→→→76→0.084632517t
20000000→1467→→→→→107246→→→→1369→→→→→98→0.0668029991
30000000→1951→→→→→152790→→→→1853→→→→→98→0.050230651cBW:&
40000000→2403→→→→→196566→→→→2300→→→→→103→0.042863088l;\`5/q
50000000→2846→→→→→239094→→→→2722→→→→→124→0.043569923{
60000000→3257→→→→→280666→→→→3126→→→→→131→0.0402210625h
70000000→3646→→→→→321468→→→→3515→→→→→131→0.035929786~#
80000000→4033→→→→→361627→→→→3892→→→→→141→0.034961567Hm9CBn
90000000→4401→→→→→401236→→→→4259→→→→→142→0.0322653948V_I
100000000→4767→→→→→440366→→→→4617→→→→→150→0.031466331nL[QTQ
110000000→5115→→→→→479074→→→→4968→→→→→147→0.028739003X<
120000000→5441→→→→→517402→→→→5312→→→→→129→0.0237088770*L
130000000→5797→→→→→555388→→→→5650→→→→→147→0.025357944Wz(
140000000→6112→→→→→593062→→→→5982→→→→→130→0.021269634y\$N0Xt
150000000→6450→→→→→630450→→→→6309→→→→→141→0.021860465Gvh,uS
160000000→6792→→→→→667574→→→→6632→→→→→160→0.023557126~0l
170000000→7113→→→→→704453→→→→6951→→→→→162→0.0227752]OX'
180000000→7446→→→→→741104→→→→7266→→→→→180→0.0241740531]4
190000000→7794→→→→→777541→→→→7577→→→→→217→0.027841936Sgq
200000000→8096→→→→→813777→→→→7884→→→→→212→0.026185771Tv2n
210000000→8400→→→→→849824→→→→8189→→→→→211→0.025119048HEt
220000000→8699→→→→→885692→→→→8490→→→→→209→0.02402575Md&q
230000000→8978→→→→→921391→→→→8789→→→→→189→0.021051459 M
240000000→9270→→→→→956928→→→→9085→→→→→185→0.01995685 h3Dh
250000000→9565→→→→→992313→→→→9379→→→→→186→0.019445896dyYRqA
260000000→9836→→→→→1027551→→→→9670→→→→→166→0.016876779FOTa
270000000→10135→→→→→1062650→→→→9959→→→→→176→0.017365565T
280000000→10431→→→→→1097615→→→→10245→→→→→186→0.017831464W,"Mt
290000000→10701→→→→→1132452→→→→10530→→→→→171→0.015979815}
300000000→10972→→→→→1167166→→→→10813→→→→→159→0.014491433Q
310000000→11280→→→→→1201762→→→→11093→→→→→187→0.016578014=j
320000000→11589→→→→→1236243→→→→11372→→→→→217→0.018724653:gJ1
330000000→11862→→→→→1270615→→→→11649→→→→→213→0.0179565C6zb
340000000→12126→→→→→1304880→→→→11925→→→→→201→0.016575952;#nw+$
350000000→12370→→→→→1339043→→→→12199→→→→→171→0.013823767lS
360000000→12632→→→→→1373107→→→→12471→→→→→161→0.012745408yg%s}
370000000→12900→→→→→1407075→→→→12742→→→→→158→0.0122480628
380000000→13164→→→→→1440950→→→→13011→→→→→153→0.011622607Zj|^Cw
390000000→13438→→→→→1474735→→→→13279→→→→→159→0.011832118l#]9A
400000000→13712→→→→→1508433→→→→13545→→→→→167→0.012179113[a3Lz
410000000→13957→→→→→1542045→→→→13810→→→→→147→0.010532349q
420000000→14247→→→→→1575575→→→→14074→→→→→173→0.012142907VWu0;
430000000→14516→→→→→1609025→→→→14337→→→→→179→0.012331221Z^|(x
440000000→14770→→→→→1642396→→→→14598→→→→→172→0.011645227j
450000000→15030→→→→→1675691→→→→14858→→→→→172→0.011443779|
460000000→15289→→→→→1708912→→→→15117→→→→→172→0.0112499181
470000000→15559→→→→→1742061→→→→15375→→→→→184→0.011825953b3:<2&
480000000→15823→→→→→1775139→→→→15632→→→→→191→0.012071036B|wI,j
490000000→16093→→→→→1808148→→→→15888→→→→→205→0.012738458u,?Pb2
500000000→16330→→→→→1841090→→→→16143→→→→→187→0.011451317A1
510000000→16579→→→→→1873967→→→→16396→→→→→183→0.01103806Sk6!kK
520000000→16816→→→→→1906779→→→→16649→→→→→167→0.009931018u}r\$b
530000000→17069→→→→→1939528→→→→16901→→→→→168→0.009842404m8Z8?