File Coverage

sieve_cluster.c
Criterion Covered Total %
statement 201 207 97.1
branch 221 288 76.7
condition n/a
subroutine n/a
pod n/a
total 422 495 85.2


line stmt bran cond sub pod time code
1             #include
2             #include
3              
4             #define FUNC_is_prime_in_sieve 1
5             #define FUNC_gcd_ui 1
6             #include "sieve.h"
7             #include "ptypes.h"
8             #include "util.h"
9             #include "primality.h"
10              
11             #define NSMALLPRIMES 168
12             static const unsigned short sprimes[NSMALLPRIMES] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941,947,953,967,971,977,983,991,997};
13              
14             typedef struct {
15             uint32_t nmax;
16             uint32_t nsize;
17             UV* list;
18             } vlist;
19             #define INIT_VLIST(v) \
20             v.nsize = 0; \
21             v.nmax = 100; \
22             New(0, v.list, v.nmax, UV);
23             #define PUSH_VLIST(v, n) \
24             do { \
25             if (v.nsize >= v.nmax) \
26             Renew(v.list, v.nmax += 100, UV); \
27             v.list[v.nsize++] = n; \
28             } while (0)
29              
30             #define ADDVAL32(v, n, max, val) \
31             do { if (n >= max) Renew(v, max += 512, UV); v[n++] = val; } while (0)
32             #define SWAPL32(l1, n1, m1, l2, n2, m2) \
33             { UV t_, *u_ = l1; l1 = l2; l2 = u_; \
34             t_ = n1; n1 = n2; n2 = t_; \
35             t_ = m1; m1 = m2; m2 = t_; }
36              
37 32           static int is_admissible(uint32_t nc, uint32_t* cl) {
38             uint32_t i, j, c;
39 32           char rset[sprimes[NSMALLPRIMES-1]];
40              
41 32 50         if (nc > NSMALLPRIMES) return 1; /* TODO */
42 241 100         for (i = 0; i < nc; i++) {
43 211           uint32_t p = sprimes[i];
44 211           memset(rset, 0, p);
45 1906 100         for (c = 0; c < nc; c++)
46 1695           rset[cl[c] % p] = 1;
47 676 100         for (j = 0; j < p; j++)
48 674 100         if (rset[j] == 0)
49 209           break;
50 211 100         if (j == p) /* All values were 1 */
51 2           return 0;
52             }
53 32           return 1;
54             }
55              
56             /* Given p prime, is this a cluster? */
57 88           static int is_cluster(UV p, uint32_t nc, uint32_t* cl) {
58             uint32_t c;
59 153 100         for (c = 1; c < nc; c++)
60 141 100         if (!is_prob_prime(p+cl[c]))
61 76           break;
62 88           return (c == nc);
63             }
64              
65             /* This is fine for small ranges. Low overhead. */
66 32           UV* sieve_cluster_simple(UV beg, UV end, uint32_t nc, uint32_t* cl, UV* numret)
67             {
68             vlist retlist;
69              
70 32 50         INIT_VLIST(retlist);
71 32 100         if (beg <= 2 && end >= 2 && is_cluster(2, nc, cl)) PUSH_VLIST(retlist, 2);
    50          
    50          
    0          
    0          
72 32 100         if (beg <= 3 && end >= 3 && is_cluster(3, nc, cl)) PUSH_VLIST(retlist, 3);
    50          
    100          
    50          
    0          
73 32 100         if (beg <= 5 && end >= 5 && is_cluster(5, nc, cl)) PUSH_VLIST(retlist, 5);
    50          
    100          
    50          
    0          
74 32 100         if (beg < 7) beg = 7;
75              
76             /* If not admissible, then don't keep looking. */
77 32 100         if (!is_admissible(nc, cl) && end > sprimes[nc])
    50          
78 2           end = sprimes[nc];
79              
80 32 50         if (beg <= end) {
81             uint32_t c;
82             unsigned char* segment;
83             UV seg_base, seg_beg, seg_end;
84              
85 32           void* ctx = start_segment_primes(beg, end, &segment);
86 67 100         while (next_segment_primes(ctx, &seg_base, &seg_beg, &seg_end)) {
87 35 100         UV sp, last_sieve_cluster = (seg_end >= cl[nc-1]) ? seg_end-cl[nc-1] : 0;
88 273643 50         START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_base, seg_beg, seg_end )
    100          
    100          
    100          
    100          
89 259643 100         if (p <= last_sieve_cluster) {
90 259591           sp = p - seg_base;
91 289914 100         for (c = 1; c < nc; c++)
92 277909 100         if (!is_prime_in_sieve(segment, sp+cl[c]))
93 247586           break;
94 259591 100         if (c == nc)
95 259591 100         PUSH_VLIST(retlist,p);
    50          
96             } else {
97 52 100         if (is_cluster(p, nc, cl))
98 2 50         PUSH_VLIST(retlist, p);
    0          
99             }
100 13963           END_DO_FOR_EACH_SIEVE_PRIME
101             }
102 32           end_segment_primes(ctx);
103             }
104 32           *numret = retlist.nsize;
105 32           return retlist.list;
106             }
107              
108              
109             #define addmodded(r,a,b,n) do { r = a + b; if (r >= n) r -= n; } while(0)
110              
111 32           UV* sieve_cluster(UV low, UV high, uint32_t nc, uint32_t* cl, UV* numret)
112             {
113             vlist retlist;
114             UV i, ppr, nres, allocres;
115 32           uint32_t const targres = 100000;
116 32           UV *residues, *cres, num_mr = 0, num_lucas = 0;
117             uint32_t pp_0, pp_1, pp_2, *resmod_0, *resmod_1, *resmod_2;
118             uint32_t rem_0, rem_1, rem_2, remadd_0, remadd_1, remadd_2;
119 32           uint32_t pi, startpi = 1, maxpi = 150;
120 32           uint32_t lastspr = sprimes[maxpi-1];
121             uint32_t c, smallnc;
122             char crem_0[43*47], crem_1[53*59], crem_2[61*67], **VPrem;
123              
124 32 50         if ((UV_MAX - cl[nc-1]) < high) return 0; /* Overflow */
125              
126 32 100         if ( ((high-low) < 10000)
127 11 100         || (nc == 3 && ((high>>31) >> 16) == 0) /* sieving large vals is slow */
    50          
128 8 100         || (nc == 2 && ((high>>31) >> 27) == 0)
    50          
129 7 50         || (nc < 2) )
130 25           return sieve_cluster_simple(low, high, nc, cl, numret);
131              
132 7 100         if (!(low&1)) low++;
133 7 50         if (!(high&1)) high--;
134              
135 7 50         INIT_VLIST(retlist);
136              
137 7 50         if (low < lastspr) {
138 7           UV t, chigh = (high > lastspr) ? lastspr : high;
139 7           UV* s = sieve_cluster_simple(low, chigh, nc, cl, &t);
140 22 100         for (i = 0; i < t; i++)
141 15 50         PUSH_VLIST(retlist, s[i]);
    0          
142 7           Safefree(s);
143 7           low = chigh + 2;
144             }
145 7 50         if (low > high) {
146 0           *numret = retlist.nsize;
147 0           return retlist.list;
148             }
149 7 50         if (low&1) low--;
150              
151             /* Determine the primorial size and acceptable residues */
152 7           New(0, residues, allocres = 1024, UV);
153             {
154             UV remr, *res2, allocres2, nres2, maxppr;
155             /* Calculate residues for a small primorial */
156 28 100         for (pi = 2, ppr = 1, i = 0; i <= pi; i++) ppr *= sprimes[i];
157 7           remr = low % ppr;
158 7           nres = 0;
159 112 100         for (i = 1; i <= ppr; i += 2) {
160 216 100         for (c = 0; c < nc; c++) {
161 210           UV v = (remr + i + cl[c]) % ppr;
162 210 100         if (gcd_ui(v, ppr) != 1) break;
163             }
164 105 100         if (c == nc)
165 6 50         ADDVAL32(residues, nres, allocres, i);
    0          
166             }
167             /* Raise primorial size until we have plenty of residues */
168 7           New(0, res2, allocres2 = 1024, UV);
169 7           maxppr = high - low;
170             #if BITS_PER_WORD == 64
171 31 50         while (pi++ < 12) {
172             #else
173             while (pi++ < 8) {
174             #endif
175 31           uint32_t j, p = sprimes[pi];
176 31           UV r, newppr = ppr * p;
177 31 100         if (nres == 0 || nres > targres/(p/2) || newppr > maxppr) break;
    50          
    100          
178 24 50         MPUverbose(2, "cluster sieve found %"UVuf" residues mod %"UVuf"\n", nres, ppr);
179 24           remr = low % newppr;
180 24           nres2 = 0;
181 312 100         for (i = 0; i < p; i++) {
182 9082 100         for (j = 0; j < nres; j++) {
183 8794           r = i*ppr + residues[j];
184 43800 100         for (c = 0; c < nc; c++) {
185 37642           UV v = remr + r + cl[c];
186 37642 100         if ((v % p) == 0) break;
187             }
188 8794 100         if (c == nc)
189 6158 100         ADDVAL32(res2, nres2, allocres2, r);
    50          
190             }
191             }
192 24           ppr = newppr;
193 24           SWAPL32(residues, nres, allocres, res2, nres2, allocres2);
194             }
195 7           startpi = pi;
196 7           Safefree(res2);
197             }
198 7 50         MPUverbose(1, "cluster sieve using %"UVuf" residues mod %"UVuf"\n", nres, ppr);
199              
200             /* Return if not admissible, maybe with a single small value */
201 7 100         if (nres == 0) {
202 1           Safefree(residues);
203 1           *numret = retlist.nsize;
204 1           return retlist.list;
205             }
206              
207             /* Pre-mod the residues with first two primes for fewer modulos every chunk */
208             {
209 6           uint32_t p1 = sprimes[startpi+0], p2 = sprimes[startpi+1];
210 6           uint32_t p3 = sprimes[startpi+2], p4 = sprimes[startpi+3];
211 6           uint32_t p5 = sprimes[startpi+4], p6 = sprimes[startpi+5];
212 6           pp_0 = p1*p2; pp_1 = p3*p4; pp_2 = p5*p6;
213 6           memset(crem_0, 1, pp_0);
214 6           memset(crem_1, 1, pp_1);
215 6           memset(crem_2, 1, pp_2);
216             /* Mark remainders that indicate a composite for this residue. */
217 120 100         for (i = 0; i < p1; i++) { crem_0[i*p1]=0; crem_0[i*p2]=0; }
218 30 100         for ( ; i < p2; i++) { crem_0[i*p1]=0; }
219 180 100         for (i = 0; i < p3; i++) { crem_1[i*p3]=0; crem_1[i*p4]=0; }
220 18 100         for ( ; i < p4; i++) { crem_1[i*p3]=0; }
221 228 100         for (i = 0; i < p5; i++) { crem_2[i*p5]=0; crem_2[i*p6]=0; }
222 30 100         for ( ; i < p6; i++) { crem_2[i*p5]=0; }
223 34 100         for (c = 1; c < nc; c++) {
224 28           uint32_t c1=cl[c], c2=cl[c], c3=cl[c], c4=cl[c], c5=cl[c], c6=cl[c];
225 28 100         if (c1 >= p1) c1 %= p1;
226 28 50         if (c2 >= p2) c2 %= p2;
227 560 100         for (i = 1; i <= p1; i++) { crem_0[i*p1-c1]=0; crem_0[i*p2-c2]=0; }
228 140 100         for ( ; i <= p2; i++) { crem_0[i*p1-c1]=0; }
229 28 50         if (c3 >= p3) c3 %= p3;
230 28 50         if (c4 >= p4) c4 %= p4;
231 840 100         for (i = 1; i <= p3; i++) { crem_1[i*p3-c3]=0; crem_1[i*p4-c4]=0; }
232 84 100         for ( ; i <= p4; i++) { crem_1[i*p3-c3]=0; }
233 28 50         if (c5 >= p5) c5 %= p5;
234 28 50         if (c6 >= p6) c6 %= p6;
235 1064 100         for (i = 1; i <= p5; i++) { crem_2[i*p5-c5]=0; crem_2[i*p6-c6]=0; }
236 140 100         for ( ; i <= p6; i++) { crem_2[i*p5-c5]=0; }
237             }
238 6 50         New(0, resmod_0, nres, uint32_t);
239 6 50         New(0, resmod_1, nres, uint32_t);
240 6 50         New(0, resmod_2, nres, uint32_t);
241 5632 100         for (i = 0; i < nres; i++) {
242 5626           resmod_0[i] = residues[i] % pp_0;
243 5626           resmod_1[i] = residues[i] % pp_1;
244 5626           resmod_2[i] = residues[i] % pp_2;
245             }
246             }
247              
248             /* Precalculate acceptable residues for more primes */
249 6 50         New(0, VPrem, maxpi, char*);
250 6           memset(VPrem, 0, maxpi);
251 828 100         for (pi = startpi+6; pi < maxpi; pi++) {
252 822           uint32_t p = sprimes[pi];
253 822           New(0, VPrem[pi], p, char);
254 822           memset(VPrem[pi], 1, p);
255             }
256 828 100         for (pi = startpi+6, smallnc = 0; pi < maxpi; pi++) {
257 822           uint32_t p = sprimes[pi];
258 822           char* prem = VPrem[pi];
259 822           prem[0] = 0;
260 856 100         while (smallnc < nc && cl[smallnc] < p) smallnc++;
    50          
261 4658 100         for (c = 1; c < smallnc; c++) prem[p-cl[c]] = 0;
262 822 50         for ( ; c < nc; c++) prem[p-(cl[c]%p)] = 0;
263             }
264              
265 6 50         New(0, cres, nres, UV);
266              
267 6           rem_0 = low % pp_0; remadd_0 = ppr % pp_0;
268 6           rem_1 = low % pp_1; remadd_1 = ppr % pp_1;
269 6           rem_2 = low % pp_2; remadd_2 = ppr % pp_2;
270              
271             /* Loop over their range in chunks of size 'ppr' */
272 24 100         while (low <= high) {
273             uint32_t r, nr, remr, ncres;
274              
275             /* Reduce the allowed residues for this chunk using more primes */
276              
277             { /* Start making a list of this chunk's residues using three pairs */
278 16896 100         for (r = 0, ncres = 0; r < nres; r++) {
279 16878 100         addmodded(remr, rem_0, resmod_0[r], pp_0);
280 16878 100         if (crem_0[remr]) {
281 9995 100         addmodded(remr, rem_1, resmod_1[r], pp_1);
282 9995 100         if (crem_1[remr]) {
283 7113 100         addmodded(remr, rem_2, resmod_2[r], pp_2);
284 7113 100         if (crem_2[remr]) {
285 5520           cres[ncres++] = residues[r];
286             }
287             }
288             }
289             }
290 18 100         addmodded(rem_0, rem_0, remadd_0, pp_0);
291 18 50         addmodded(rem_1, rem_1, remadd_1, pp_1);
292 18 100         addmodded(rem_2, rem_2, remadd_2, pp_2);
293             }
294              
295             /* Sieve through more primes one at a time, removing residues. */
296 2455 100         for (pi = startpi+6; pi < maxpi && ncres > 0; pi++) {
    100          
297 2437           uint32_t p = sprimes[pi];
298 2437           uint32_t rem = low % p;
299 2437           char* prem = VPrem[pi];
300             /* Check divisibility of each remaining residue with this p */
301             /* If we extended prem we could remove the add in the loop below */
302 2437 50         if (startpi <= 9) { /* Residues are 32-bit */
303 142234 100         for (r = 0, nr = 0; r < ncres; r++) {
304 139797 100         if (prem[ (rem+(uint32_t)cres[r]) % p ])
305 134664           cres[nr++] = cres[r];
306             }
307             } else { /* Residues are 64-bit */
308 0 0         for (r = 0, nr = 0; r < ncres; r++) {
309 0 0         if (prem[ (rem+cres[r]) % p ])
310 0           cres[nr++] = cres[r];
311             }
312             }
313 2437           ncres = nr;
314             }
315 18 50         MPUverbose(3, "cluster sieve range has %u residues left\n", ncres);
316              
317             /* Now check each of the remaining residues for inclusion */
318 286 100         for (r = 0; r < ncres; r++) {
319 272           UV p = low + cres[r];
320 272 100         if (p > high) break;
321             /* PRP test. Split to save time. */
322 1425 100         for (c = 0; c < nc; c++)
323 1166 100         if (num_mr++,!is_euler_plumb_pseudoprime(p+cl[c]))
324 9           break;
325 268 100         if (c < nc) continue;
326 1393 100         for (c = 0; c < nc; c++)
327 1134 50         if (num_lucas++,!is_almost_extra_strong_lucas_pseudoprime(p+cl[c], 1))
328 0           break;
329 259 50         if (c < nc) continue;
330 259 100         PUSH_VLIST(retlist, p);
    50          
331             }
332 18           low += ppr;
333 18 50         if (low < ppr) low = UV_MAX;
334             }
335              
336 6 50         MPUverbose(1, "cluster sieve ran %"UVuf" MR and %"UVuf" Lucas tests\n", num_mr, num_lucas);
337 828 100         for (pi = startpi+6; pi < maxpi; pi++)
338 822           Safefree(VPrem[pi]);
339 6           Safefree(VPrem);
340 6           Safefree(resmod_0);
341 6           Safefree(resmod_1);
342 6           Safefree(resmod_2);
343 6           Safefree(cres);
344 6           Safefree(residues);
345 6           *numret = retlist.nsize;
346 32           return retlist.list;
347             }