File Coverage

sieve.c
Criterion Covered Total %
statement 181 200 90.5
branch 269 340 79.1
condition n/a
subroutine n/a
pod n/a
total 450 540 83.3


line stmt bran cond sub pod time code
1             #include
2             #include
3             #include
4             #include
5              
6             #define FUNC_isqrt 1
7             #define FUNC_next_prime_in_sieve
8             #include "sieve.h"
9             #include "ptypes.h"
10             #include "cache.h"
11             #include "util.h"
12             #include "primality.h"
13             #include "montmath.h"
14             #include "prime_nth_count.h"
15              
16             /* Is it better to do a partial sieve + primality tests vs. full sieve? */
17 1870           static int do_partial_sieve(UV startp, UV endp) {
18 1870           UV range = endp - startp;
19 1870           if (USE_MONTMATH) range /= 8; /* Fast primality tests */
20             #if BITS_PER_WORD == 64
21 1870 100         if ( (startp > UVCONST( 100000000000000) && range < 40000) ||
    50          
    50          
22 1859 0         (startp > UVCONST( 1000000000000000) && range < 150000) ||
    50          
23 1859 0         (startp > UVCONST( 10000000000000000) && range < 600000) ||
    50          
24 1859 0         (startp > UVCONST( 100000000000000000) && range < 2500000) ||
    50          
25 1859 0         (startp > UVCONST( 1000000000000000000) && range < 10000000) ||
    50          
26 0 0         (startp > UVCONST(10000000000000000000) && range < 40000000) )
27 11           return 1;
28             #endif
29 1859           return 0;
30             }
31              
32             /* 1001 bytes of presieved mod-30 bytes. If the area to be sieved is
33             * appropriately filled with this data, then 7, 11, and 13 do not have
34             * to be sieved. It wraps, so multiple memcpy's can be used. Do be
35             * aware that if you start at 0, you'll have to correct the first byte.
36             *
37             * mpu '$g=7*11*13; @b=(0)x$g; for $d (0..$g-1) { $i=0; for $m (1,7,11,13,17,19,23,29) { $n=30*$d+$m; if (gcd($n,$g) != 1) { $b[$d] |= (1<<$i); } $i++; } } for (0..$#b) { printf "0x%02x,",$b[$_]; print "\n" unless ($_+1)%13; } print "\n"'
38             */
39             #define PRESIEVE_SIZE (7*11*13)
40             static const unsigned char presieve13[PRESIEVE_SIZE] =
41             { 0x0e,0x20,0x10,0x81,0x49,0x24,0xc2,0x06,0x2a,0x90,0xa1,0x0c,0x14,
42             0x58,0x02,0x61,0x11,0xc3,0x28,0x0c,0x44,0x22,0xa4,0x10,0x91,0x18,
43             0x4d,0x40,0x82,0x21,0x58,0xa1,0x28,0x04,0x42,0x92,0x20,0x51,0x91,
44             0x8a,0x04,0x48,0x03,0x60,0x34,0x81,0x1c,0x06,0xc1,0x02,0xa2,0x10,
45             0x89,0x08,0x24,0x45,0x42,0x30,0x10,0xc5,0x0a,0x86,0x40,0x0a,0x30,
46             0x38,0x85,0x08,0x15,0x40,0x63,0x20,0x96,0x83,0x88,0x04,0x60,0x16,
47             0x28,0x10,0x81,0x49,0x44,0xe2,0x02,0x2c,0x12,0xa1,0x0c,0x04,0x50,
48             0x0a,0x61,0x10,0x83,0x48,0x2c,0x40,0x26,0x26,0x90,0x91,0x08,0x55,
49             0x48,0x82,0x20,0x19,0xc1,0x28,0x04,0x44,0x12,0xa0,0x51,0x81,0x9a,
50             0x0c,0x48,0x02,0x21,0x54,0xa1,0x18,0x04,0x43,0x82,0xa2,0x10,0x99,
51             0x08,0x24,0x44,0x03,0x70,0x30,0xc1,0x0c,0x86,0xc0,0x0a,0x20,0x30,
52             0x8d,0x08,0x14,0x41,0x43,0x20,0x92,0x85,0x0a,0x84,0x60,0x06,0x30,
53             0x18,0x81,0x49,0x05,0xc2,0x22,0x28,0x14,0xa3,0x8c,0x04,0x50,0x12,
54             0x69,0x10,0x83,0x09,0x4c,0x60,0x22,0x24,0x12,0x91,0x08,0x45,0x50,
55             0x8a,0x20,0x18,0x81,0x68,0x24,0x40,0x16,0x22,0xd1,0x81,0x8a,0x14,
56             0x48,0x02,0x20,0x15,0xc1,0x38,0x04,0x45,0x02,0xa2,0x10,0x89,0x18,
57             0x2c,0x44,0x02,0x31,0x50,0xe1,0x08,0x86,0x42,0x8a,0x20,0x30,0x95,
58             0x08,0x14,0x40,0x43,0x60,0xb2,0x81,0x0c,0x06,0xe0,0x06,0x20,0x10,
59             0x89,0x49,0x04,0xc3,0x42,0x28,0x10,0xa5,0x0e,0x84,0x50,0x02,0x71,
60             0x18,0x83,0x08,0x0d,0x40,0x22,0x24,0x14,0x93,0x88,0x45,0x40,0x92,
61             0x28,0x18,0x81,0x29,0x44,0x60,0x12,0x24,0x53,0x81,0x8a,0x04,0x58,
62             0x0a,0x20,0x14,0x81,0x58,0x24,0x41,0x06,0xa2,0x90,0x89,0x08,0x34,
63             0x4c,0x02,0x30,0x11,0xc1,0x28,0x86,0x44,0x0a,0xa0,0x30,0x85,0x18,
64             0x1c,0x40,0x43,0x21,0xd2,0xa1,0x08,0x04,0x62,0x86,0x20,0x10,0x91,
65             0x49,0x04,0xc2,0x03,0x68,0x30,0xa1,0x0c,0x06,0xd0,0x02,0x61,0x10,
66             0x8b,0x08,0x0c,0x41,0x62,0x24,0x10,0x95,0x0a,0xc5,0x40,0x82,0x30,
67             0x18,0x81,0x28,0x05,0x40,0x32,0x20,0x55,0x83,0x8a,0x04,0x48,0x12,
68             0x28,0x14,0x81,0x19,0x44,0x61,0x02,0xa6,0x12,0x89,0x08,0x24,0x54,
69             0x0a,0x30,0x10,0xc1,0x48,0xa6,0x40,0x0e,0x22,0xb0,0x85,0x08,0x14,
70             0x48,0x43,0x20,0x93,0xc1,0x28,0x04,0x64,0x06,0xa0,0x10,0x81,0x59,
71             0x0c,0xc2,0x02,0x29,0x50,0xa1,0x0c,0x04,0x52,0x82,0x61,0x10,0x93,
72             0x08,0x0c,0x40,0x23,0x64,0x30,0x91,0x0c,0x47,0xc0,0x82,0x20,0x18,
73             0x89,0x28,0x04,0x41,0x52,0x20,0x51,0x85,0x8a,0x84,0x48,0x02,0x30,
74             0x1c,0x81,0x18,0x05,0x41,0x22,0xa2,0x14,0x8b,0x88,0x24,0x44,0x12,
75             0x38,0x10,0xc1,0x09,0xc6,0x60,0x0a,0x24,0x32,0x85,0x08,0x14,0x50,
76             0x4b,0x20,0x92,0x81,0x48,0x24,0x60,0x06,0x22,0x90,0x81,0x49,0x14,
77             0xca,0x02,0x28,0x11,0xe1,0x2c,0x04,0x54,0x02,0xe1,0x10,0x83,0x18,
78             0x0c,0x40,0x22,0x25,0x50,0xb1,0x08,0x45,0x42,0x82,0x20,0x18,0x91,
79             0x28,0x04,0x40,0x13,0x60,0x71,0x81,0x8e,0x06,0xc8,0x02,0x20,0x14,
80             0x89,0x18,0x04,0x41,0x42,0xa2,0x10,0x8d,0x0a,0xa4,0x44,0x02,0x30,
81             0x18,0xc1,0x08,0x87,0x40,0x2a,0x20,0x34,0x87,0x88,0x14,0x40,0x53,
82             0x28,0x92,0x81,0x09,0x44,0x60,0x06,0x24,0x12,0x81,0x49,0x04,0xd2,
83             0x0a,0x28,0x10,0xa1,0x4c,0x24,0x50,0x06,0x63,0x90,0x83,0x08,0x1c,
84             0x48,0x22,0x24,0x11,0xd1,0x28,0x45,0x44,0x82,0xa0,0x18,0x81,0x38,
85             0x0c,0x40,0x12,0x21,0x51,0xa1,0x8a,0x04,0x4a,0x82,0x20,0x14,0x91,
86             0x18,0x04,0x41,0x03,0xe2,0x30,0x89,0x0c,0x26,0xc4,0x02,0x30,0x10,
87             0xc9,0x08,0x86,0x41,0x4a,0x20,0x30,0x85,0x0a,0x94,0x40,0x43,0x30,
88             0x9a,0x81,0x08,0x05,0x60,0x26,0x20,0x14,0x83,0xc9,0x04,0xc2,0x12,
89             0x28,0x10,0xa1,0x0d,0x44,0x70,0x02,0x65,0x12,0x83,0x08,0x0c,0x50,
90             0x2a,0x24,0x10,0x91,0x48,0x65,0x40,0x86,0x22,0x98,0x81,0x28,0x14,
91             0x48,0x12,0x20,0x51,0xc1,0xaa,0x04,0x4c,0x02,0xa0,0x14,0x81,0x18,
92             0x0c,0x41,0x02,0xa3,0x50,0xa9,0x08,0x24,0x46,0x82,0x30,0x10,0xd1,
93             0x08,0x86,0x40,0x0b,0x60,0x30,0x85,0x0c,0x16,0xc0,0x43,0x20,0x92,
94             0x89,0x08,0x04,0x61,0x46,0x20,0x10,0x85,0x4b,0x84,0xc2,0x02,0x38,
95             0x18,0xa1,0x0c,0x05,0x50,0x22,0x61,0x14,0x83,0x88,0x0c,0x40,0x32,
96             0x2c,0x10,0x91,0x09,0x45,0x60,0x82,0x24,0x1a,0x81,0x28,0x04,0x50,
97             0x1a,0x20,0x51,0x81,0xca,0x24,0x48,0x06,0x22,0x94,0x81,0x18,0x14,
98             0x49,0x02,0xa2,0x11,0xc9,0x28,0x24,0x44,0x02,0xb0,0x10,0xc1,0x18,
99             0x8e,0x40,0x0a,0x21,0x70,0xa5,0x08,0x14,0x42,0xc3,0x20,0x92,0x91,
100             0x08,0x04,0x60,0x07,0x60,0x30,0x81,0x4d,0x06,0xc2,0x02,0x28,0x10,
101             0xa9,0x0c,0x04,0x51,0x42,0x61,0x10,0x87,0x0a,0x8c,0x40,0x22,0x34,
102             0x18,0x91,0x08,0x45,0x40,0xa2,0x20,0x1c,0x83,0xa8,0x04,0x40,0x12,
103             0x28,0x51,0x81,0x8b,0x44,0x68,0x02,0x24,0x16,0x81,0x18,0x04,0x51,
104             0x0a,0xa2,0x10,0x89,0x48,0x24,0x44,0x06,0x32,0x90,0xc1,0x08,0x96,
105             0x48,0x0a,0x20,0x31,0xc5,0x28,0x14,0x44,0x43,0xa0,0x92,0x81,0x18,
106             0x0c,0x60,0x06,0x21,0x50,0xa1,0x49,0x04,0xc2,0x82,0x28,0x10,0xb1,
107             0x0c,0x04,0x50,0x03,0x61,0x30,0x83,0x0c,0x0e,0xc0,0x22,0x24,0x10,
108             0x99,0x08,0x45,0x41,0xc2,0x20,0x18,0x85,0x2a,0x84,0x40,0x12,0x30,
109             0x59,0x81,0x8a,0x05,0x48,0x22,0x20,0x14,0x83,0x98,0x04,0x41,0x12,
110             0xaa,0x10,0x89,0x09,0x64,0x64,0x02,0x34,0x12,0xc1,0x08,0x86,0x50,
111             0x0a,0x20,0x30,0x85,0x48,0x34,0x40,0x47,0x22,0x92,0x81,0x08,0x14,
112             0x68,0x06,0x20,0x11,0xc1,0x69,0x04,0xc6,0x02,0xa8,0x10,0xa1,0x1c,
113             0x0c,0x50,0x02,0x61,0x50,0xa3,0x08,0x0c,0x42,0xa2,0x24,0x10,0x91,
114             0x08,0x45,0x40,0x83,0x60,0x38,0x81,0x2c,0x06,0xc0,0x12,0x20,0x51,
115             0x89,0x8a,0x04,0x49,0x42,0x20,0x14,0x85,0x1a,0x84,0x41,0x02,0xb2,
116             0x18,0x89,0x08,0x25,0x44,0x22,0x30,0x14,0xc3,0x88,0x86,0x40,0x1a,
117             0x28,0x30,0x85,0x09,0x54,0x60,0x43,0x24,0x92,0x81,0x08,0x04,0x70};
118              
119             static const UV max_sieve_prime = (BITS_PER_WORD==64) ? 4294967291U : 65521U;
120              
121             /* Tile bytes from source to bytes in dest */
122 1251           static void memtile(unsigned char* dst, const unsigned char* src, size_t from, size_t to) {
123 1251 100         if (to < from)
124 1092           from = to;
125 1251 50         if (dst != src)
126 0           memcpy(dst, src, from);
127 1898 100         while (from < to) {
128 647 100         size_t bytes = (2*from > to) ? to-from : from;
129 647           memcpy(dst+from, dst, bytes);
130 647           from += bytes;
131             }
132 1251           }
133              
134 2016           static UV sieve_prefill(unsigned char* mem, UV startd, UV endd)
135             {
136 2016           UV vnext_prime = 17;
137 2016           UV nbytes = endd - startd + 1;
138 2016 50         MPUassert( (mem != 0) && (endd >= startd), "sieve_prefill bad arguments");
    50          
139              
140 2016 100         if (startd != 0) {
141 819           UV pstartd = startd % PRESIEVE_SIZE;
142 819           UV tailbytes = PRESIEVE_SIZE - pstartd;
143 819 100         if (tailbytes > nbytes) tailbytes = nbytes;
144 819           memcpy(mem, presieve13 + pstartd, tailbytes); /* Copy tail to mem */
145 819           mem += tailbytes; /* Advance so mem points at the beginning */
146 819           nbytes -= tailbytes;
147             }
148 2016 100         if (nbytes > 0) {
149 1251           memcpy(mem, presieve13, (nbytes < PRESIEVE_SIZE) ? nbytes : PRESIEVE_SIZE);
150 1251           memtile(mem, mem, PRESIEVE_SIZE, nbytes);
151 1251 100         if (startd == 0) mem[0] = 0x01; /* Correct first byte */
152             }
153             /* Leaving option open to tile 17 out and sieve, then return 19 */
154 2016           return vnext_prime;
155             }
156              
157             /* Marking primes is done the same way we used to do with tables, but
158             * now uses heavily unrolled code based on Kim Walisch's mod-30 sieve.
159             */
160             #define set_bit(s,n) *(s) |= (1 << n);
161             static const unsigned char masknum30[30] =
162             {0,0,0,0,0,0,0,1,0,0,0,2,0,3,0,0,0,4,0,5,0,0,0,6,0,0,0,0,0,7};
163             static const unsigned char qinit30[30] =
164             {0,0,1,1,1,1,1,1,2,2,2,2,3,3,4,4,4,4,5,5,6,6,6,6,7,7,7,7,7,7};
165              
166             typedef struct {
167             uint32_t prime;
168             uint32_t offset;
169             uint8_t index;
170             } wheel_t;
171              
172             #define CROSS_INDEX(v, b0,b1,b2,b3,b4,b5,b6,b7, i0,i1,i2,i3,i4,i5,i6,i7, it) \
173             while (1) { \
174             case (v+0): if(s>=send){w->index=v+0;break;} set_bit(s,b0); s += r*6+i0; \
175             case (v+1): if(s>=send){w->index=v+1;break;} set_bit(s,b1); s += r*4+i1; \
176             case (v+2): if(s>=send){w->index=v+2;break;} set_bit(s,b2); s += r*2+i2; \
177             case (v+3): if(s>=send){w->index=v+3;break;} set_bit(s,b3); s += r*4+i3; \
178             case (v+4): if(s>=send){w->index=v+4;break;} set_bit(s,b4); s += r*2+i4; \
179             case (v+5): if(s>=send){w->index=v+5;break;} set_bit(s,b5); s += r*4+i5; \
180             case (v+6): if(s>=send){w->index=v+6;break;} set_bit(s,b6); s += r*6+i6; \
181             case (v+7): if(s>=send){w->index=v+7;break;} set_bit(s,b7); s += r*2+i7; \
182             while (s + r*28 + it-1 < send) { \
183             set_bit(s + r * 0 + 0, b0); \
184             set_bit(s + r * 6 + i0, b1); \
185             set_bit(s + r * 10 + i0+i1, b2); \
186             set_bit(s + r * 12 + i0+i1+i2, b3); \
187             set_bit(s + r * 16 + i0+i1+i2+i3, b4); \
188             set_bit(s + r * 18 + i0+i1+i2+i3+i4, b5); \
189             set_bit(s + r * 22 + i0+i1+i2+i3+i4+i5, b6); \
190             set_bit(s + r * 28 + i0+i1+i2+i3+i4+i5+i6, b7); \
191             s += r*30 + it; \
192             } \
193             }
194 4236785           static wheel_t create_wheel(UV startp, uint32_t prime)
195             {
196             wheel_t w;
197 4236785           UV q = prime;
198 4236785           UV p2 = q*q;
199              
200 4236785 100         if (startp == 0) {
201 15881           wheel_t ws = { prime, p2/30, qinit30[q % 30] + 8*masknum30[prime % 30] };
202 15881           return ws;
203             }
204              
205 4220904 100         if (p2 < startp) {
206 4218061           q = 1+(startp-1)/prime;
207 4218061           q += distancewheel30[q % 30];
208 4218061           p2 = prime * q;
209             /* The offset if p2 overflows is still ok, or set to max_sieve_prime+1. */
210             /* if (p2 < startp) p2 = max_sieve_prime+1; */
211             }
212              
213 4220904           w.offset = (p2-startp) / 30;
214 4220904           w.index = qinit30[q % 30] + 8*masknum30[prime % 30];
215 4220904           w.prime = prime;
216 4236785           return w;
217             }
218              
219 4243940           static void mark_primes(unsigned char* s, uint32_t bytes, wheel_t* w)
220             {
221 4243940 100         if (w->offset >= bytes) {
222 4098422           w->offset -= bytes;
223             } else {
224 145518           const unsigned char* send = s + bytes;
225 145518           uint32_t r = w->prime / 30;
226 145518           s += w->offset;
227 145518           switch (w->index) {
228 548192 100         CROSS_INDEX( 0, 0,1,2,3,4,5,6,7, 0,0,0,0,0,0,0,1, 1); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
229 577458 100         CROSS_INDEX( 8, 1,5,4,0,7,3,2,6, 1,1,1,0,1,1,1,1, 7); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
230 548155 100         CROSS_INDEX(16, 2,4,0,6,1,7,3,5, 2,2,0,2,0,2,2,1, 11); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
231 528945 100         CROSS_INDEX(24, 3,0,6,5,2,1,7,4, 3,1,1,2,1,1,3,1, 13); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
232 796934 100         CROSS_INDEX(32, 4,7,1,2,5,6,0,3, 3,3,1,2,1,3,3,1, 17); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
233 661452 100         CROSS_INDEX(40, 5,3,7,1,6,0,4,2, 4,2,2,2,2,2,4,1, 19); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
234 711409 100         CROSS_INDEX(48, 6,2,3,7,0,4,5,1, 5,3,1,4,1,3,5,1, 23); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
235 625933 100         CROSS_INDEX(56, 7,6,5,4,3,2,1,0, 6,4,2,4,2,4,6,1, 29); break;
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    100          
236             }
237 145518           w->offset = s - send;
238             }
239 4243940           }
240              
241             /* Monolithic mod-30 wheel sieve */
242 119           unsigned char* sieve_erat30(UV end)
243             {
244             unsigned char *mem;
245             UV max_buf, limit, prime;
246              
247 119           max_buf = (end/30) + ((end%30) != 0);
248             /* Round up to a word */
249 119           max_buf = ((max_buf + sizeof(UV) - 1) / sizeof(UV)) * sizeof(UV);
250 119           New(0, mem, max_buf, unsigned char );
251              
252             /* Fill buffer with marked 7, 11, and 13 */
253 119           prime = sieve_prefill(mem, 0, max_buf-1);
254              
255 119           limit = isqrt(end); /* prime*prime can overflow */
256 11865 100         for ( ; prime <= limit; prime = next_prime_in_sieve(mem,prime,end)) {
257 11746           wheel_t w = create_wheel(0, prime);
258 11746           mark_primes(mem, max_buf, &w);
259             }
260 119           return mem;
261             }
262              
263             /* Segmented mod-30 wheel sieve */
264 269           int sieve_segment(unsigned char* mem, UV startd, UV endd)
265             {
266             const unsigned char* sieve;
267             UV limit, slimit, start_base_prime, sieve_size;
268 269           UV startp = 30*startd;
269 269 50         UV endp = (endd >= (UV_MAX/30)) ? UV_MAX-2 : 30*endd+29;
270 269 50         MPUassert( (mem != 0) && (endd >= startd) && (endp >= startp),
    50          
    50          
271             "sieve_segment bad arguments");
272              
273             /* It's possible we can just use the primary cache */
274 269           sieve_size = get_prime_cache(0, &sieve);
275 269 100         if (sieve_size >= endp) {
276 246           memcpy(mem, sieve+startd, endd-startd+1);
277             release_prime_cache(sieve);
278 246           return 1;
279             }
280              
281             /* Fill buffer with marked 7, 11, and 13 */
282 23           start_base_prime = sieve_prefill(mem, startd, endd);
283              
284 23           limit = isqrt(endp);
285 23 50         if (limit > max_sieve_prime) limit = max_sieve_prime;
286 23           slimit = limit;
287 23 50         if (do_partial_sieve(startp, endp))
288 0 0         slimit >>= ((startp < (UV)1e16) ? 8 : 10);
289             /* printf("segment sieve from %"UVuf" to %"UVuf" (aux sieve to %"UVuf")\n", startp, endp, slimit); */
290 23 50         if (slimit > sieve_size) {
291             release_prime_cache(sieve);
292 0           get_prime_cache(slimit, &sieve);
293             }
294              
295 80936 50         START_DO_FOR_EACH_SIEVE_PRIME(sieve, 0, start_base_prime, slimit) {
    100          
    100          
    100          
    100          
296 77136           wheel_t w = create_wheel(startp, p);
297 77136           mark_primes(mem, endd-startd+1, &w);
298 3708           } END_DO_FOR_EACH_SIEVE_PRIME;
299             release_prime_cache(sieve);
300              
301 23 50         if (limit > slimit) { /* We've sieved out most composites, but not all. */
302 0 0         START_DO_FOR_EACH_SIEVE_PRIME(mem, 0, 0, endp-startp) {
    0          
    0          
    0          
    0          
303 0 0         if (!BPSW(startp + p)) /* If the candidate is not prime, */
304 0           mem[p/30] |= masktab30[p%30]; /* mark the sieve location. */
305 0           } END_DO_FOR_EACH_SIEVE_PRIME;
306             }
307 269           return 1;
308             }
309              
310 1874           int sieve_segment_wheel(unsigned char* mem, UV startd, UV endd, wheel_t *warray, uint32_t wsize)
311             {
312 1874           uint32_t i = 0, limit, start_base_prime;
313 1874           uint32_t segsize = endd - startd + 1;
314 1874           UV startp = 30*startd;
315 1874 50         UV endp = (endd >= (UV_MAX/30)) ? UV_MAX-2 : 30*endd+29;
316 1874 50         MPUassert( (mem != 0) && (endd >= startd) && (endp >= startp),
    50          
    50          
317             "sieve_segment bad arguments");
318              
319             /* possibly use primary cache directly */
320              
321             /* Fill buffer with marked 7, 11, and 13 */
322 1874           start_base_prime = sieve_prefill(mem, startd, endd);
323 12884 100         while (i < wsize && warray[i].prime < start_base_prime)
    100          
324 11010           i++;
325              
326 1874           limit = isqrt(endp);
327 1874 50         if (limit > max_sieve_prime) limit = max_sieve_prime;
328              
329 4156932 100         while (i < wsize && warray[i].prime <= limit) {
    100          
330 4155058 100         if (warray[i].index >= 64)
331 4147903           warray[i] = create_wheel(startp, warray[i].prime);
332 4155058           mark_primes(mem, segsize, &(warray[i++]));
333             }
334              
335 1874 100         if (limit > warray[wsize-1].prime && warray[wsize-1].prime < max_sieve_prime) {
    50          
336 7487 50         START_DO_FOR_EACH_SIEVE_PRIME(mem, 0, 0, endp-startp) {
    100          
    50          
    100          
    100          
337 6871 100         if (!BPSW(startp + p)) /* If the candidate is not prime, */
338 2789           mem[p/30] |= masktab30[p%30]; /* mark the sieve location. */
339 605           } END_DO_FOR_EACH_SIEVE_PRIME;
340             }
341              
342 1874           return 1;
343             }
344              
345             /**************************************************************************/
346              
347             typedef struct {
348             UV lod;
349             UV hid;
350             UV low;
351             UV high;
352             UV endp;
353             UV segment_size;
354             unsigned char* segment;
355             unsigned char* base;
356             wheel_t *warray;
357             uint32_t wsize;
358             } segment_context_t;
359              
360             /*
361             * unsigned char* segment;
362             * UV seg_base, seg_low, seg_high;
363             * void* ctx = start_segment_primes(low, high, &segment);
364             * while (beg < 7) {
365             * beg = (beg <= 2) ? 2 : (beg <= 3) ? 3 : 5;
366             * .... with beg ....
367             * beg += 1 + (beg > 2);
368             * }
369             * while (next_segment_primes(ctx, &seg_base, &seg_low, &seg_high)) {
370             * START_DO_FOR_EACH_SIEVE_PRIME( segment, seg_low - seg_base, seg_high - seg_base )
371             * .... with seg_base + p ....
372             * END_DO_FOR_EACH_SIEVE_PRIME
373             * }
374             * end_segment_primes(ctx);
375             */
376              
377 1847           void* start_segment_primes(UV low, UV high, unsigned char** segmentmem)
378             {
379             segment_context_t* ctx;
380             UV nsegments, range;
381              
382 1847 50         MPUassert( high >= low, "start_segment_primes bad arguments");
383 1847           New(0, ctx, 1, segment_context_t);
384 1847           ctx->low = low;
385 1847           ctx->high = high;
386 1847           ctx->lod = low / 30;
387 1847           ctx->hid = high / 30;
388 1847 50         ctx->endp = (ctx->hid >= (UV_MAX/30)) ? UV_MAX-2 : 30*ctx->hid+29;
389 1847           range = ctx->hid - ctx->lod + 1; /* range in bytes */
390              
391             #if BITS_PER_WORD == 64
392 1847 100         if (high > 1e10 && range > 32*1024-16) {
    50          
393             UV size, div;
394             /* Use larger segments */
395 0           size = isqrt(32*isqrt(high)) * (logint(high,2)-2);
396 0 0         if (size < 128*1024) size = 128*1024;
397             /* Evenly split the range into segments */
398 0           div = (range+size-1)/size;
399 0 0         size = (div <= 1) ? range : (range+div-1)/div;
400 0           ctx->segment_size = size;
401 0           New(0, ctx->segment, size, unsigned char);
402             } else
403             #endif
404 1847           ctx->segment = get_prime_segment( &(ctx->segment_size) );
405 1847           *segmentmem = ctx->segment;
406 1847           nsegments = (((high-low+29)/30)+ctx->segment_size-1) / ctx->segment_size;
407              
408 1847 50         if (_XS_get_verbose() >= 2)
409 0           printf("segment sieve: byte range %lu split into %lu segments of size %lu\n", (unsigned long)range, (unsigned long)nsegments, (unsigned long)ctx->segment_size);
410              
411 1847           ctx->base = 0;
412 1847           ctx->warray = 0;
413 1847           ctx->wsize = 0;
414             #if 1
415             { /* Generate wheel data for this segment sieve */
416 1847           const UV maxsieve = UVCONST(400000000);
417             UV limit, nprimes;
418             wheel_t *warray;
419 1847           wheel_t w = {0,0,128};
420 1847           uint32_t wsize = 0;
421             /* Number of primes for a full sieve */
422 1847           limit = isqrt(ctx->endp);
423             /* For small ranges a partial sieve is much faster */
424 1847 100         if (do_partial_sieve(low, high))
425 11 100         limit >>= ((low < (UV)1e16) ? 8 : 10);
426 1847 50         if (limit <= maxsieve) {
427             /* Bump to one more than needed. */
428 1847           limit = next_prime(limit);
429             /* We'll make space for this many */
430 1847           nprimes = prime_count_upper(limit);
431 1847 50         if (_XS_get_verbose() >= 4)
432 0           printf("segment sieve %lu - %lu, primes to %lu (max %lu)\n", (unsigned long)low, (unsigned long)high, (unsigned long)limit, (unsigned long)nprimes);
433 1847 50         New(0, warray, nprimes, wheel_t);
434 4317082 50         START_DO_FOR_EACH_PRIME(0,limit) {
    100          
    100          
    100          
    100          
    100          
    100          
    100          
    50          
    100          
435 4160348 50         if (wsize >= nprimes) croak("segment bad upper count");
436 4160348           w.prime = p;
437 4160348           warray[wsize++] = w;
438 4160348           } END_DO_FOR_EACH_PRIME;
439 1847           ctx->warray = warray;
440 1847           ctx->wsize = wsize;
441             }
442             }
443             #endif
444              
445 1847           return (void*) ctx;
446             }
447              
448 3648           int next_segment_primes(void* vctx, UV* base, UV* low, UV* high)
449             {
450             UV seghigh_d, range_d;
451 3648           segment_context_t* ctx = (segment_context_t*) vctx;
452              
453 3648 100         if (ctx->lod > ctx->hid) return 0;
454              
455 3748           seghigh_d = ((ctx->hid - ctx->lod) < ctx->segment_size)
456             ? ctx->hid
457 1874 100         : (ctx->lod + ctx->segment_size - 1);
458 1874           range_d = seghigh_d - ctx->lod + 1;
459 1874           *low = ctx->low;
460 1874 100         *high = (seghigh_d == ctx->hid) ? ctx->high : (seghigh_d*30 + 29);
461 1874           *base = ctx->lod * 30;
462              
463 1874 50         MPUassert( seghigh_d >= ctx->lod, "next_segment_primes: highd < lowd");
464 1874 50         MPUassert( range_d <= ctx->segment_size, "next_segment_primes: range > segment size");
465              
466 1874 50         if (ctx->warray != 0)
467 1874           sieve_segment_wheel(ctx->segment, ctx->lod, seghigh_d, ctx->warray, ctx->wsize);
468             else
469 0           sieve_segment(ctx->segment, ctx->lod, seghigh_d);
470              
471 1874           ctx->lod += range_d;
472 1874           ctx->low = *high + 2;
473              
474 1874           return 1;
475             }
476              
477 1847           void end_segment_primes(void* vctx)
478             {
479 1847           segment_context_t* ctx = (segment_context_t*) vctx;
480 1847 50         MPUassert(ctx != 0, "end_segment_primes given a null pointer");
481 1847 50         if (ctx->segment != 0) {
482 1847           release_prime_segment(ctx->segment);
483 1847           ctx->segment = 0;
484             }
485 1847 50         if (ctx->base != 0) {
486 0           Safefree(ctx->base);
487 0           ctx->base = 0;
488             }
489 1847 50         if (ctx->warray != 0) {
490 1847           Safefree(ctx->warray);
491 1847           ctx->warray = 0;
492             }
493 1847           Safefree(ctx);
494 1847           }