1 /*
2 Copyright 2018-2019 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 [n0] <nMax>
31
32 Checks mpz_nextprime() exhaustively, starting from n=n0 up to
33 nMax.
34
35 WARNING: The full intervall [0..nMax] is sieved at once, even if
36 only a piece is needed. This may require a lot of memory!
37
38 */
39
40 #include <stdlib.h>
41 #include <stdio.h>
42 #include "gmp-impl.h"
43 #include "longlong.h"
44 #include "tests.h"
45 #define STOP(x) return (x)
46 /* #define STOP(x) x */
47 #define REPS 10
48 /* #define TRACE(x,n) if ((n)>1) {x;} */
49 #define TRACE(x,n)
50
51 /* The full primesieve.c is included, just for block_resieve, that
52 is not exported ... */
53 #undef gmp_primesieve
54 #include "../../primesieve.c"
55
56 #ifndef BLOCK_SIZE
57 #define BLOCK_SIZE 2048
58 #endif
59
60 /*********************************************************/
61 /* Section sieve: sieving functions and tools for primes */
62 /*********************************************************/
63
64 static mp_size_t
primesieve_size(mp_limb_t n)65 primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; }
66
67 /*************************************************************/
68 /* Section macros: common macros, for swing/fac/bin (&sieve) */
69 /*************************************************************/
70
71 #define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve) \
72 __max_i = (end); \
73 \
74 do { \
75 ++__i; \
76 if (((sieve)[__index] & __mask) == 0) \
77 { \
78 mp_limb_t prime; \
79 prime = id_to_n(__i)
80
81 #define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve) \
82 do { \
83 mp_limb_t __mask, __index, __max_i, __i; \
84 \
85 __i = (start)-(off); \
86 __index = __i / GMP_LIMB_BITS; \
87 __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS); \
88 __i += (off); \
89 \
90 LOOP_ON_SIEVE_CONTINUE(prime,end,sieve)
91
92 #define LOOP_ON_SIEVE_STOP \
93 } \
94 __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1); \
95 __index += __mask & 1; \
96 } while (__i <= __max_i)
97
98 #define LOOP_ON_SIEVE_END \
99 LOOP_ON_SIEVE_STOP; \
100 } while (0)
101
102 mpz_t g;
103
something_wrong(mpz_t er,int exp)104 int something_wrong (mpz_t er, int exp)
105 {
106 fprintf (stderr, "value = %lu , expected = %i\n", mpz_get_ui (er), exp);
107 return -1;
108 }
109
110 int
check_pprime(unsigned long begin,unsigned long end,int composites)111 check_pprime (unsigned long begin, unsigned long end, int composites)
112 {
113 begin = (begin / 6U) * 6U;
114 for (;(begin < 2) & (begin <= end); ++begin)
115 {
116 *(g->_mp_d) = begin;
117 TRACE(printf ("-%li ", begin),1);
118 if (mpz_probab_prime_p (g, REPS))
119 STOP (something_wrong (g, 0));
120 }
121 for (;(begin < 4) & (begin <= end); ++begin)
122 {
123 *(g->_mp_d) = begin;
124 TRACE(printf ("+%li ", begin),2);
125 if (!composites && !mpz_probab_prime_p (g, REPS))
126 STOP (something_wrong (g, 1));
127 }
128 if (end > 4) {
129 if ((end > 10000) && (begin > end / begin))
130 {
131 mp_limb_t *sieve, *primes;
132 mp_size_t size_s, size_p, off;
133 unsigned long start;
134
135 mpz_set_ui (g, end);
136 mpz_sqrt (g, g);
137 start = mpz_get_ui (g) + GMP_LIMB_BITS;
138 size_p = primesieve_size (start);
139
140 primes = __GMP_ALLOCATE_FUNC_LIMBS (size_p);
141 gmp_primesieve (primes, start);
142
143 size_s = BLOCK_SIZE * 2;
144 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size_s);
145 off = n_to_bit(begin) + (begin % 3 == 0);
146
147 do {
148 TRACE (printf ("off =%li\n", off),3);
149 block_resieve (sieve, BLOCK_SIZE, off, primes);
150 TRACE (printf ("LOOP =%li - %li\n", id_to_n (off+1), id_to_n (off + BLOCK_SIZE * GMP_LIMB_BITS)),3);
151 LOOP_ON_SIEVE_BEGIN (prime, off, off + BLOCK_SIZE * GMP_LIMB_BITS - 1,
152 off, sieve);
153
154 do {
155 *(g->_mp_d) = begin;
156 TRACE(printf ("-%li ", begin),1);
157 if (mpz_probab_prime_p (g, REPS))
158 STOP (something_wrong (g, 0));
159 if ((begin & 0xff) == 0)
160 {
161 spinner();
162 if ((begin & 0xfffffff) == 0)
163 printf ("%li (0x%lx)\n", begin, begin);
164 }
165 } while (++begin < prime);
166
167 *(g->_mp_d) = begin;
168 TRACE(printf ("+%li ", begin),2);
169 if (!composites && ! mpz_probab_prime_p (g, REPS))
170 STOP (something_wrong (g, 1));
171 ++begin;
172
173 LOOP_ON_SIEVE_END;
174 off += BLOCK_SIZE * GMP_LIMB_BITS;
175 } while (begin < end);
176
177 __GMP_FREE_FUNC_LIMBS (sieve, size_s);
178 __GMP_FREE_FUNC_LIMBS (primes, size_p);
179 }
180 else
181 {
182 mp_limb_t *sieve;
183 mp_size_t size;
184 unsigned long start;
185
186 size = primesieve_size (end);
187
188 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
189 gmp_primesieve (sieve, end);
190 start = MAX (begin, 5) | 1;
191 LOOP_ON_SIEVE_BEGIN (prime, n_to_bit(start) + (start % 3 == 0),
192 n_to_bit (end), 0, sieve);
193
194 do {
195 *(g->_mp_d) = begin;
196 TRACE(printf ("-%li ", begin),1);
197 if (mpz_probab_prime_p (g, REPS))
198 STOP (something_wrong (g, 0));
199 if ((begin & 0xff) == 0)
200 {
201 spinner();
202 if ((begin & 0xfffffff) == 0)
203 printf ("%li (0x%lx)\n", begin, begin);
204 }
205 } while (++begin < prime);
206
207 *(g->_mp_d) = begin;
208 TRACE(printf ("+%li ", begin),2);
209 if (!composites && ! mpz_probab_prime_p (g, REPS))
210 STOP (something_wrong (g, 1));
211 ++begin;
212
213 LOOP_ON_SIEVE_END;
214
215 __GMP_FREE_FUNC_LIMBS (sieve, size);
216 }
217 }
218
219 for (;begin < end; ++begin)
220 {
221 *(g->_mp_d) = begin;
222 TRACE(printf ("-%li ", begin),1);
223 if (mpz_probab_prime_p (g, REPS))
224 STOP (something_wrong (g, 0));
225 }
226
227 gmp_printf ("%Zd\n", g);
228 return 0;
229 }
230
231 int
check_nprime(unsigned long begin,unsigned long end)232 check_nprime (unsigned long begin, unsigned long end)
233 {
234 if (begin < 2)
235 {
236 *(g->_mp_d) = begin;
237 TRACE(printf ("%li ", begin),1);
238 mpz_nextprime (g, g);
239 if (mpz_cmp_ui (g, 2) != 0)
240 STOP (something_wrong (g, -1));
241 begin = mpz_get_ui (g);
242 }
243 if (begin < 3)
244 {
245 *(g->_mp_d) = begin;
246 TRACE(printf ("%li ", begin),1);
247 mpz_nextprime (g, g);
248 if (mpz_cmp_ui (g, 3) != 0)
249 STOP (something_wrong (g, -1));
250 begin = mpz_get_ui (g);
251 }
252 if (end > 4)
253 {
254 mp_limb_t *sieve;
255 mp_size_t size;
256 unsigned long start;
257
258 size = primesieve_size (end);
259
260 sieve = __GMP_ALLOCATE_FUNC_LIMBS (size);
261 gmp_primesieve (sieve, end);
262 start = MAX (begin, 5) | 1;
263 *(g->_mp_d) = begin;
264 LOOP_ON_SIEVE_BEGIN (prime, n_to_bit(start) + (start % 3 == 0),
265 n_to_bit (end), 0, sieve);
266
267 mpz_nextprime (g, g);
268 if (mpz_cmp_ui (g, prime) != 0)
269 STOP (something_wrong (g, -1));
270
271 if (prime - start > 200)
272 {
273 start = prime;
274 spinner();
275 if (prime - begin > 0xfffffff)
276 {
277 begin = prime;
278 printf ("%li (0x%lx)\n", begin, begin);
279 }
280 }
281
282 LOOP_ON_SIEVE_END;
283
284 __GMP_FREE_FUNC_LIMBS (sieve, size);
285 }
286
287 if (mpz_cmp_ui (g, end) < 0)
288 {
289 mpz_nextprime (g, g);
290 if (mpz_cmp_ui (g, end) <= 0)
291 STOP (something_wrong (g, -1));
292 }
293
294 gmp_printf ("%Zd\n", g);
295 return 0;
296 }
297
298 int
main(int argc,char ** argv)299 main (int argc, char **argv)
300 {
301 int ret, mode = 0;
302 unsigned long begin = 0, end = 0;
303
304 for (;argc > 1;--argc,++argv)
305 switch (*argv[1]) {
306 case 'p':
307 mode = 0;
308 break;
309 case 'c':
310 mode = 2;
311 break;
312 case 'n':
313 mode = 1;
314 break;
315 default:
316 begin = end;
317 end = atol (argv[1]);
318 }
319
320 if (begin >= end)
321 {
322 fprintf (stderr, "usage: primes [n|p|c] [n0] <nMax>\n");
323 exit (1);
324 }
325
326 mpz_init_set_ui (g, ULONG_MAX);
327
328 switch (mode) {
329 case 1:
330 ret = check_nprime (begin, end);
331 break;
332 default:
333 ret = check_pprime (begin, end, mode);
334 }
335
336 mpz_clear (g);
337
338 if (ret == 0)
339 printf ("Prime tests checked in [%lu - %lu] [0x%lx - 0x%lx].\n", begin, end, begin, end);
340 return ret;
341 }
342