1 /*
2 Copyright 2018-2020 Free Software Foundation, Inc.
3
4 This file is part of the GNU MP Library test suite.
5
6 The GNU MP Library test suite is free software; you can redistribute it
7 and/or modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 3 of the License,
9 or (at your option) any later version.
10
11 The GNU MP Library test suite is distributed in the hope that it will be
12 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
14 Public License for more details.
15
16 You should have received a copy of the GNU General Public License along with
17 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
18
19 /* Usage:
20
21 ./primes [p|c] [n0] <nMax>
22
23 Checks mpz_probab_prime_p(n, r) exhaustively, starting from n=n0
24 up to nMax.
25 If n0 * n0 > nMax, the intervall is sieved piecewise, else the
26 full intervall [0..nMax] is sieved at once.
27 With the parameter "p" (or nothing), tests all numbers. With "c"
28 only composites are tested.
29
30 ./primes n|N [n0] <nMax>
31
32 Checks mpz_nextprime() exhaustively, starting from n=n0 up to
33 nMax. With "n", only the sequence of primes is checked, with "N"
34 the function is tested on every number in the interval.
35
36 WARNING: The full intervall [0..nMax] is sieved at once, even if
37 only a piece is needed. This may require a lot of memory!
38
39 */
40
41 #include <stdlib.h>
42 #include <stdio.h>
43 #include "gmp-impl.h"
44 #include "longlong.h"
45 #include "tests.h"
46 #define STOP(x) return (x)
47 /* #define STOP(x) x */
48 #define REPS 10
49 /* #define TRACE(x,n) if ((n)>1) {x;} */
50 #define TRACE(x,n)
51
52 /* The full primesieve.c is included, just for block_resieve, that
53 is not exported ... */
54 #undef gmp_primesieve
55 #include "../../primesieve.c"
56
57 #ifndef BLOCK_SIZE
58 #define BLOCK_SIZE 2048
59 #endif
60
61 /*********************************************************/
62 /* Section sieve: sieving functions and tools for primes */
63 /*********************************************************/
64
65 static mp_size_t
66 primesieve_size (mp_limb_t n) { return n_fto_bit(n) / GMP_LIMB_BITS + 1; }
67
68 /*************************************************************/
69 /* Section macros: common macros, for swing/fac/bin (&sieve) */
70 /*************************************************************/
71
72 #define LOOP_ON_SIEVE_CONTINUE(prime,end) \
73 __max_i = (end); \
74 \
75 do { \
76 ++__i; \
77 if ((*__sieve & __mask) == 0) \
78 { \
79 mp_limb_t prime; \
80 prime = id_to_n(__i)
81
82 #define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve) \
83 do { \
84 mp_limb_t __mask, *__sieve, __max_i, __i; \
85 \
86 __i = (start)-(off); \
87 __sieve = (sieve) + __i / GMP_LIMB_BITS; \
88 __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS); \
89 __i += (off); \
90 \
91 LOOP_ON_SIEVE_CONTINUE(prime,end)
92
93 #define LOOP_ON_SIEVE_STOP \
94 } \
95 __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1); \
96 __sieve += __mask & 1; \
97 } while (__i <= __max_i)
98
99 #define LOOP_ON_SIEVE_END \
100 LOOP_ON_SIEVE_STOP; \
101 } while (0)
102
103 mpz_t g;
104
105 int something_wrong (mpz_t er, int exp)
106 {
107 fprintf (stderr, "value = %lu , expected = %i\n", mpz_get_ui (er), exp);
108 return -1;
109 }
110
111 int
112 check_pprime (unsigned long begin, unsigned long end, int composites)
113 {
114 begin = (begin / 6U) * 6U;
115 for (;(begin < 2) & (begin <= end); ++begin)
116 {
117 *(g->_mp_d) = begin;
118 TRACE(printf ("-%li ", begin),1);
119 if (mpz_probab_prime_p (g, REPS))
120 STOP (something_wrong (g, 0));
121 }
122 for (;(begin < 4) & (begin <= end); ++begin)
123 {
124 *(g->_mp_d) = begin;
125 TRACE(printf ("+%li ", begin),2);
126 if (!composites && !mpz_probab_prime_p (g, REPS))
127 STOP (something_wrong (g, 1));
128 }
129 if (end > 4) {
130 if ((end > 10000) && (begin > end / begin))
131 {
132 mp_limb_t *sieve, *primes;
133 mp_size_t size_s, size_p, off;
134 unsigned long start;
135
136 mpz_set_ui (g, end);
137 mpz_sqrt (g, g);
138 start = mpz_get_ui (g) + GMP_LIMB_BITS;
139 size_p = primesieve_size (start);
140
141 primes = __GMP_ALLOCATE_FUNC_LIMBS (size_p);
142 gmp_primesieve (primes, start);
143
144 size_s = BLOCK_SIZE * 2;
145 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size_s);
146 off = n_cto_bit(begin);
147
148 do {
149 TRACE (printf ("off =%li\n", off),3);
150 block_resieve (sieve, BLOCK_SIZE, off, primes);
151 TRACE (printf ("LOOP =%li - %li\n", id_to_n (off+1), id_to_n (off + BLOCK_SIZE * GMP_LIMB_BITS)),3);
152 LOOP_ON_SIEVE_BEGIN (prime, off, off + BLOCK_SIZE * GMP_LIMB_BITS - 1,
153 off, sieve);
154
155 do {
156 *(g->_mp_d) = begin;
157 TRACE(printf ("-%li ", begin),1);
158 if (mpz_probab_prime_p (g, REPS))
159 STOP (something_wrong (g, 0));
160 if ((begin & 0xff) == 0)
161 {
162 spinner();
163 if ((begin & 0xfffffff) == 0)
164 printf ("%li (0x%lx)\n", begin, begin);
165 }
166 } while (++begin < prime);
167
168 *(g->_mp_d) = begin;
169 TRACE(printf ("+%li ", begin),2);
170 if (!composites && ! mpz_probab_prime_p (g, REPS))
171 STOP (something_wrong (g, 1));
172 ++begin;
173
174 LOOP_ON_SIEVE_END;
175 off += BLOCK_SIZE * GMP_LIMB_BITS;
176 } while (begin < end);
177
178 __GMP_FREE_FUNC_LIMBS (sieve, size_s);
179 __GMP_FREE_FUNC_LIMBS (primes, size_p);
180 }
181 else
182 {
183 mp_limb_t *sieve;
184 mp_size_t size;
185 unsigned long start;
186
187 size = primesieve_size (end);
188
189 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
190 gmp_primesieve (sieve, end);
191 start = MAX (begin, 5) | 1;
192 LOOP_ON_SIEVE_BEGIN (prime, n_cto_bit(start),
193 n_fto_bit (end), 0, sieve);
194
195 do {
196 *(g->_mp_d) = begin;
197 TRACE(printf ("-%li ", begin),1);
198 if (mpz_probab_prime_p (g, REPS))
199 STOP (something_wrong (g, 0));
200 if ((begin & 0xff) == 0)
201 {
202 spinner();
203 if ((begin & 0xfffffff) == 0)
204 printf ("%li (0x%lx)\n", begin, begin);
205 }
206 } while (++begin < prime);
207
208 *(g->_mp_d) = begin;
209 TRACE(printf ("+%li ", begin),2);
210 if (!composites && ! mpz_probab_prime_p (g, REPS))
211 STOP (something_wrong (g, 1));
212 ++begin;
213
214 LOOP_ON_SIEVE_END;
215
216 __GMP_FREE_FUNC_LIMBS (sieve, size);
217 }
218 }
219
220 for (;begin < end; ++begin)
221 {
222 *(g->_mp_d) = begin;
223 TRACE(printf ("-%li ", begin),1);
224 if (mpz_probab_prime_p (g, REPS))
225 STOP (something_wrong (g, 0));
226 }
227
228 gmp_printf ("%Zd\n", g);
229 return 0;
230 }
231
232 int
233 check_nprime (unsigned long begin, unsigned long end)
234 {
235 if (begin < 2)
236 {
237 *(g->_mp_d) = begin;
238 g->_mp_size = begin;
239 TRACE(printf ("%li ", begin),1);
240 mpz_nextprime (g, g);
241 if (mpz_cmp_ui (g, 2) != 0)
242 STOP (something_wrong (g, 2));
243 begin = mpz_get_ui (g);
244 }
245 if (begin < 3)
246 {
247 *(g->_mp_d) = begin;
248 TRACE(printf ("%li ", begin),1);
249 mpz_nextprime (g, g);
250 if (mpz_cmp_ui (g, 3) != 0)
251 STOP (something_wrong (g, 3));
252 begin = mpz_get_ui (g);
253 }
254 if (end > 4)
255 {
256 mp_limb_t *sieve;
257 mp_size_t size;
258 unsigned long start;
259
260 size = primesieve_size (end);
261
262 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
263 gmp_primesieve (sieve, end);
264 start = MAX (begin, 5) | 1;
265 *(g->_mp_d) = begin;
266 LOOP_ON_SIEVE_BEGIN (prime, n_cto_bit(start),
267 n_fto_bit (end), 0, sieve);
268
269 mpz_nextprime (g, g);
270 if (mpz_cmp_ui (g, prime) != 0)
271 STOP (something_wrong (g, prime));
272
273 if (prime - start > 200)
274 {
275 start = prime;
276 spinner();
277 if (prime - begin > 0xfffffff)
278 {
279 begin = prime;
280 printf ("%li (0x%lx)\n", begin, begin);
281 }
282 }
283
284 LOOP_ON_SIEVE_END;
285
286 __GMP_FREE_FUNC_LIMBS (sieve, size);
287 }
288
289 if (mpz_cmp_ui (g, end) < 0)
290 {
291 mpz_nextprime (g, g);
292 if (mpz_cmp_ui (g, end) <= 0)
293 STOP (something_wrong (g, -1));
294 }
295
296 gmp_printf ("%Zd\n", g);
297 return 0;
298 }
299
300 int
301 check_Nprime (unsigned long begin, unsigned long end)
302 {
303 mpz_t op;
304 mpz_init_set_ui (op, end);
305
306 for (;begin < 2; ++begin)
307 {
308 *(op->_mp_d) = begin;
309 op->_mp_size = begin;
310 TRACE(printf ("%li ", begin),1);
311 mpz_nextprime (g, op);
312 if (mpz_cmp_ui (g, 2) != 0)
313 STOP (something_wrong (g, 2));
314 }
315 if (begin < 3)
316 {
317 *(op->_mp_d) = begin;
318 TRACE(printf ("%li ", begin),1);
319 mpz_nextprime (g, op);
320 if (mpz_cmp_ui (g, 3) != 0)
321 STOP (something_wrong (g, 3));
322 begin = 3;
323 }
324 if (end > 4)
325 {
326 mp_limb_t *sieve;
327 mp_size_t size;
328 unsigned long start;
329 unsigned long opl;
330
331 size = primesieve_size (end);
332
333 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
334 gmp_primesieve (sieve, end);
335 start = MAX (begin, 5) | 1;
336 opl = begin;
337 LOOP_ON_SIEVE_BEGIN (prime, n_cto_bit(start),
338 n_fto_bit (end), 0, sieve);
339
340 do {
341 *(op->_mp_d) = opl;
342 mpz_nextprime (g, op);
343 if (mpz_cmp_ui (g, prime) != 0)
344 STOP (something_wrong (g, prime));
345 ++opl;
346 } while (opl < prime);
347
348 if (prime - start > 200)
349 {
350 start = prime;
351 spinner();
352 if (prime - begin > 0xfffffff)
353 {
354 begin = prime;
355 printf ("%li (0x%lx)\n", begin, begin);
356 }
357 }
358
359 LOOP_ON_SIEVE_END;
360
361 __GMP_FREE_FUNC_LIMBS (sieve, size);
362 }
363
364 if (mpz_cmp_ui (g, end) < 0)
365 {
366 mpz_nextprime (g, g);
367 if (mpz_cmp_ui (g, end) <= 0)
368 STOP (something_wrong (g, -1));
369 }
370
371 gmp_printf ("%Zd\n", g);
372 return 0;
373 }
374
375 int
376 main (int argc, char **argv)
377 {
378 int ret, mode = 0;
379 unsigned long begin = 0, end = 0;
380
381 for (;argc > 1;--argc,++argv)
382 switch (*argv[1]) {
383 case 'p':
384 mode = 0;
385 break;
386 case 'c':
387 mode = 2;
388 break;
389 case 'n':
390 mode = 1;
391 break;
392 case 'N':
393 mode = 3;
394 break;
395 default:
396 begin = end;
397 end = atol (argv[1]);
398 }
399
400 if (begin >= end)
401 {
402 fprintf (stderr, "usage: primes [N|n|p|c] [n0] <nMax>\n");
403 exit (1);
404 }
405
406 mpz_init_set_ui (g, ULONG_MAX);
407
408 switch (mode) {
409 case 1:
410 ret = check_nprime (begin, end);
411 break;
412 case 3:
413 ret = check_Nprime (begin, end);
414 break;
415 default:
416 ret = check_pprime (begin, end, mode);
417 }
418
419 mpz_clear (g);
420
421 if (ret == 0)
422 printf ("Prime tests checked in [%lu - %lu] [0x%lx - 0x%lx].\n", begin, end, begin, end);
423 return ret;
424 }