1 /* primesieve (BIT_ARRAY, N) -- Fills the BIT_ARRAY with a mask for primes up to N.
2
3 Contributed to the GNU project by Marco Bodrato.
4
5 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.
6 IT IS ONLY SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.
7 IN FACT, IT IS ALMOST GUARANTEED THAT IT WILL CHANGE OR
8 DISAPPEAR IN A FUTURE GNU MP RELEASE.
9
10 Copyright 2010-2012, 2015, 2016, 2021, 2022 Free Software Foundation, Inc.
11
12 This file is part of the GNU MP Library.
13
14 The GNU MP Library is free software; you can redistribute it and/or modify
15 it under the terms of either:
16
17 * the GNU Lesser General Public License as published by the Free
18 Software Foundation; either version 3 of the License, or (at your
19 option) any later version.
20
21 or
22
23 * the GNU General Public License as published by the Free Software
24 Foundation; either version 2 of the License, or (at your option) any
25 later version.
26
27 or both in parallel, as here.
28
29 The GNU MP Library is distributed in the hope that it will be useful, but
30 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
31 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
32 for more details.
33
34 You should have received copies of the GNU General Public License and the
35 GNU Lesser General Public License along with the GNU MP Library. If not,
36 see https://www.gnu.org/licenses/. */
37
38 #include "gmp-impl.h"
39
40 #if 0
41 static mp_limb_t
42 bit_to_n (mp_limb_t bit) { return (bit*3+4)|1; }
43 #endif
44
45 /* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/
46 static mp_limb_t
47 id_to_n (mp_limb_t id) { return id*3+1+(id&1); }
48
49 /* n_fto_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */
50 static mp_limb_t
51 n_fto_bit (mp_limb_t n) { return ((n-5)|1)/3U; }
52
53 /* n_cto_bit (n) = ((n-2)&(-CNST_LIMB(2)))/3U */
54 static mp_limb_t
55 n_cto_bit (mp_limb_t n) { return (n|1)/3U-1; }
56
57 #if 0
58 static mp_size_t
59 primesieve_size (mp_limb_t n) { return n_fto_bit(n) / GMP_LIMB_BITS + 1; }
60 #endif
61
62 #define SET_OFF1(m1, m2, M1, M2, off, BITS) \
63 if (off) { \
64 if (off < GMP_LIMB_BITS) { \
65 m1 = (M1 >> off) | (M2 << (GMP_LIMB_BITS - off)); \
66 if (off <= BITS - GMP_LIMB_BITS) { \
67 m2 = M1 << (BITS - GMP_LIMB_BITS - off) \
68 | M2 >> off; \
69 } else { \
70 m1 |= M1 << (BITS - off); \
71 m2 = M1 >> (off + GMP_LIMB_BITS - BITS); \
72 } \
73 } else { \
74 m1 = M1 << (BITS - off) \
75 | M2 >> (off - GMP_LIMB_BITS); \
76 m2 = M2 << (BITS - off) \
77 | M1 >> (off + GMP_LIMB_BITS - BITS); \
78 } \
79 } else { \
80 m1 = M1; m2 = M2; \
81 }
82
83 #define SET_OFF2(m1, m2, m3, M1, M2, M3, off, BITS) \
84 if (off) { \
85 if (off <= GMP_LIMB_BITS) { \
86 m1 = M2 << (GMP_LIMB_BITS - off); \
87 m2 = M3 << (GMP_LIMB_BITS - off); \
88 if (off != GMP_LIMB_BITS) { \
89 m1 |= (M1 >> off); \
90 m2 |= (M2 >> off); \
91 } \
92 if (off <= BITS - 2 * GMP_LIMB_BITS) { \
93 m3 = M1 << (BITS - 2 * GMP_LIMB_BITS - off) \
94 | M3 >> off; \
95 } else { \
96 m2 |= M1 << (BITS - GMP_LIMB_BITS - off); \
97 m3 = M1 >> (off + 2 * GMP_LIMB_BITS - BITS); \
98 } \
99 } else if (off < 2 *GMP_LIMB_BITS) { \
100 m1 = M2 >> (off - GMP_LIMB_BITS) \
101 | M3 << (2 * GMP_LIMB_BITS - off); \
102 if (off <= BITS - GMP_LIMB_BITS) { \
103 m2 = M3 >> (off - GMP_LIMB_BITS) \
104 | M1 << (BITS - GMP_LIMB_BITS - off); \
105 m3 = M2 << (BITS - GMP_LIMB_BITS - off); \
106 if (off != BITS - GMP_LIMB_BITS) { \
107 m3 |= M1 >> (off + 2 * GMP_LIMB_BITS - BITS); \
108 } \
109 } else { \
110 m1 |= M1 << (BITS - off); \
111 m2 = M2 << (BITS - off) \
112 | M1 >> (GMP_LIMB_BITS - BITS + off); \
113 m3 = M2 >> (GMP_LIMB_BITS - BITS + off); \
114 } \
115 } else { \
116 m1 = M1 << (BITS - off) \
117 | M3 >> (off - 2 * GMP_LIMB_BITS); \
118 m2 = M2 << (BITS - off) \
119 | M1 >> (off + GMP_LIMB_BITS - BITS); \
120 m3 = M3 << (BITS - off) \
121 | M2 >> (off + GMP_LIMB_BITS - BITS); \
122 } \
123 } else { \
124 m1 = M1; m2 = M2; m3 = M3; \
125 }
126
127 #define ROTATE1(m1, m2, BITS) \
128 do { \
129 mp_limb_t __tmp; \
130 __tmp = m1 >> (2 * GMP_LIMB_BITS - BITS); \
131 m1 = (m1 << (BITS - GMP_LIMB_BITS)) | m2; \
132 m2 = __tmp; \
133 } while (0)
134
135 #define ROTATE2(m1, m2, m3, BITS) \
136 do { \
137 mp_limb_t __tmp; \
138 __tmp = m2 >> (3 * GMP_LIMB_BITS - BITS); \
139 m2 = m2 << (BITS - GMP_LIMB_BITS * 2) \
140 | m1 >> (3 * GMP_LIMB_BITS - BITS); \
141 m1 = m1 << (BITS - GMP_LIMB_BITS * 2) | m3; \
142 m3 = __tmp; \
143 } while (0)
144
145 static mp_limb_t
146 fill_bitpattern (mp_ptr bit_array, mp_size_t limbs, mp_limb_t offset)
147 {
148 #ifdef SIEVE_2MSK2
149 mp_limb_t m11, m12, m21, m22, m23;
150
151 { /* correctly handle offset == 0... */
152 mp_limb_t off1 = offset % (11 * 5 * 2);
153 SET_OFF1 (m11, m12, SIEVE_MASK1, SIEVE_MASKT, off1, 11 * 5 * 2);
154 offset %= 13 * 7 * 2;
155 SET_OFF2 (m21, m22, m23, SIEVE_2MSK1, SIEVE_2MSK2, SIEVE_2MSKT, offset, 13 * 7 * 2);
156 }
157 /* THINK: Consider handling odd values of 'limbs' outside the loop,
158 to have a single exit condition. */
159 do {
160 bit_array[0] = m11 | m21;
161 if (--limbs == 0)
162 break;
163 ROTATE1 (m11, m12, 11 * 5 * 2);
164 bit_array[1] = m11 | m22;
165 bit_array += 2;
166 ROTATE1 (m11, m12, 11 * 5 * 2);
167 ROTATE2 (m21, m22, m23, 13 * 7 * 2);
168 } while (--limbs != 0);
169 return n_cto_bit (13 + 1);
170 #else
171 #ifdef SIEVE_MASK2
172 mp_limb_t mask, mask2, tail;
173
174 { /* correctly handle offset == 0... */
175 offset %= 7 * 5 * 2;
176 SET_OFF2 (mask, mask2, tail, SIEVE_MASK1, SIEVE_MASK2, SIEVE_MASKT, offset, 7 * 5 * 2);
177 }
178 /* THINK: Consider handling odd values of 'limbs' outside the loop,
179 to have a single exit condition. */
180 do {
181 bit_array[0] = mask;
182 if (--limbs == 0)
183 break;
184 bit_array[1] = mask2;
185 bit_array += 2;
186 ROTATE2 (mask, mask2, tail, 7 * 5 * 2);
187 } while (--limbs != 0);
188 return n_cto_bit (7 + 1);
189 #else
190 MPN_FILL (bit_array, limbs, CNST_LIMB(0));
191 return 0;
192 #endif
193 #endif
194 }
195
196 static void
197 block_resieve (mp_ptr bit_array, mp_size_t limbs, mp_limb_t offset,
198 mp_srcptr sieve)
199 {
200 mp_size_t bits, off = offset;
201 mp_limb_t mask, i;
202
203 ASSERT (limbs > 0);
204
205 bits = limbs * GMP_LIMB_BITS - 1;
206
207 i = fill_bitpattern (bit_array, limbs, offset);
208
209 ASSERT (i < GMP_LIMB_BITS);
210
211 mask = CNST_LIMB(1) << i;
212 do {
213 ++i;
214 if ((*sieve & mask) == 0)
215 {
216 mp_size_t step, lindex;
217 mp_limb_t lmask;
218 unsigned maskrot;
219
220 step = id_to_n(i);
221
222 /* lindex = n_to_bit(id_to_n(i)*id_to_n(i)); */
223 lindex = i*(step+1)-1+(-(i&1)&(i+1));
224 /* lindex = i*(step+1+(i&1))-1+(i&1); */
225 if (lindex > bits + off)
226 break;
227
228 step <<= 1;
229 maskrot = step % GMP_LIMB_BITS;
230
231 if (lindex < off)
232 lindex += step * ((off - lindex - 1) / step + 1);
233
234 lindex -= off;
235
236 lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
237 for ( ; lindex <= bits; lindex += step) {
238 bit_array[lindex / GMP_LIMB_BITS] |= lmask;
239 lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
240 };
241
242 /* lindex = n_to_bit(id_to_n(i)*bit_to_n(i)); */
243 lindex = i*(i*3+6)+(i&1);
244
245 if (lindex < off)
246 lindex += step * ((off - lindex - 1) / step + 1);
247
248 lindex -= off;
249
250 lmask = CNST_LIMB(1) << (lindex % GMP_LIMB_BITS);
251 for ( ; lindex <= bits; lindex += step) {
252 bit_array[lindex / GMP_LIMB_BITS] |= lmask;
253 lmask = lmask << maskrot | lmask >> (GMP_LIMB_BITS - maskrot);
254 };
255 }
256 mask = mask << 1 | mask >> (GMP_LIMB_BITS-1);
257 sieve += mask & 1;
258 } while (1);
259 }
260
261 #define BLOCK_SIZE 2048
262
263 /* Fills bit_array with the characteristic function of composite
264 numbers up to the parameter n. I.e. a bit set to "1" represents a
265 composite, a "0" represents a prime.
266
267 The primesieve_size(n) limbs pointed to by bit_array are
268 overwritten. The returned value counts prime integers in the
269 interval [4, n]. Note that n > 4.
270
271 Even numbers and multiples of 3 are excluded "a priori", only
272 numbers equivalent to +/- 1 mod 6 have their bit in the array.
273
274 Once sieved, if the bit b is ZERO it represent a prime, the
275 represented prime is bit_to_n(b), if the LSbit is bit 0, or
276 id_to_n(b), if you call "1" the first bit.
277 */
278
279 mp_limb_t
280 gmp_primesieve (mp_ptr bit_array, mp_limb_t n)
281 {
282 mp_size_t size;
283 mp_limb_t bits;
284 static mp_limb_t presieved[] = {PRIMESIEVE_INIT_TABLE};
285
286 ASSERT (n > 4);
287
288 bits = n_fto_bit(n);
289 size = bits / GMP_LIMB_BITS + 1;
290
291 for (mp_size_t j = 0, lim = MIN (size, PRIMESIEVE_NUMBEROF_TABLE);
292 j < lim; ++j)
293 bit_array [j] = presieved [j]; /* memcopy? */
294
295 if (size > PRIMESIEVE_NUMBEROF_TABLE) {
296 mp_size_t off;
297 off = size > 2 * BLOCK_SIZE ? BLOCK_SIZE + (size % BLOCK_SIZE) : size;
298 block_resieve (bit_array + PRIMESIEVE_NUMBEROF_TABLE,
299 off - PRIMESIEVE_NUMBEROF_TABLE,
300 GMP_LIMB_BITS * PRIMESIEVE_NUMBEROF_TABLE, bit_array);
301 for (; off < size; off += BLOCK_SIZE)
302 block_resieve (bit_array + off, BLOCK_SIZE, off * GMP_LIMB_BITS, bit_array);
303 }
304
305 if ((bits + 1) % GMP_LIMB_BITS != 0)
306 bit_array[size-1] |= MP_LIMB_T_MAX << ((bits + 1) % GMP_LIMB_BITS);
307
308 return size * GMP_LIMB_BITS - mpn_popcount (bit_array, size);
309 }
310
311 #undef BLOCK_SIZE
312 #undef SIEVE_MASK1
313 #undef SIEVE_MASK2
314 #undef SIEVE_MASKT
315 #undef SIEVE_2MSK1
316 #undef SIEVE_2MSK2
317 #undef SIEVE_2MSKT
318 #undef SET_OFF1
319 #undef SET_OFF2
320 #undef ROTATE1
321 #undef ROTATE2