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 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 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 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 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 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