1/* mpz_nextprime(p,t) - compute the next prime > t and store that in p. 2 3Copyright 1999, 2000, 2001, 2008, 2009 Free Software Foundation, Inc. 4 5Contributed to the GNU project by Niels M�ller and Torbjorn Granlund. 6 7This file is part of the GNU MP Library. 8 9The GNU MP Library is free software; you can redistribute it and/or modify 10it under the terms of the GNU Lesser General Public License as published by 11the Free Software Foundation; either version 3 of the License, or (at your 12option) any later version. 13 14The GNU MP Library is distributed in the hope that it will be useful, but 15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 16or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 17License for more details. 18 19You should have received a copy of the GNU Lesser General Public License 20along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 21 22#include "gmp.h" 23#include "gmp-impl.h" 24#include "longlong.h" 25 26static const unsigned char primegap[] = 27{ 28 2,2,4,2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,14,4,6, 29 2,10,2,6,6,4,6,6,2,10,2,4,2,12,12,4,2,4,6,2,10,6,6,6,2,6,4,2,10,14,4,2, 30 4,14,6,10,2,4,6,8,6,6,4,6,8,4,8,10,2,10,2,6,4,6,8,4,2,4,12,8,4,8,4,6, 31 12,2,18,6,10,6,6,2,6,10,6,6,2,6,6,4,2,12,10,2,4,6,6,2,12,4,6,8,10,8,10,8, 32 6,6,4,8,6,4,8,4,14,10,12,2,10,2,4,2,10,14,4,2,4,14,4,2,4,20,4,8,10,8,4,6, 33 6,14,4,6,6,8,6,12 34}; 35 36#define NUMBER_OF_PRIMES 167 37 38void 39mpz_nextprime (mpz_ptr p, mpz_srcptr n) 40{ 41 unsigned short *moduli; 42 unsigned long difference; 43 int i; 44 unsigned prime_limit; 45 unsigned long prime; 46 int cnt; 47 mp_size_t pn; 48 mp_bitcnt_t nbits; 49 unsigned incr; 50 TMP_SDECL; 51 52 /* First handle tiny numbers */ 53 if (mpz_cmp_ui (n, 2) < 0) 54 { 55 mpz_set_ui (p, 2); 56 return; 57 } 58 mpz_add_ui (p, n, 1); 59 mpz_setbit (p, 0); 60 61 if (mpz_cmp_ui (p, 7) <= 0) 62 return; 63 64 pn = SIZ(p); 65 count_leading_zeros (cnt, PTR(p)[pn - 1]); 66 nbits = pn * GMP_NUMB_BITS - (cnt - GMP_NAIL_BITS); 67 if (nbits / 2 >= NUMBER_OF_PRIMES) 68 prime_limit = NUMBER_OF_PRIMES - 1; 69 else 70 prime_limit = nbits / 2; 71 72 TMP_SMARK; 73 74 /* Compute residues modulo small odd primes */ 75 moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short); 76 77 for (;;) 78 { 79 /* FIXME: Compute lazily? */ 80 prime = 3; 81 for (i = 0; i < prime_limit; i++) 82 { 83 moduli[i] = mpz_fdiv_ui (p, prime); 84 prime += primegap[i]; 85 } 86 87#define INCR_LIMIT 0x10000 /* deep science */ 88 89 for (difference = incr = 0; incr < INCR_LIMIT; difference += 2) 90 { 91 /* First check residues */ 92 prime = 3; 93 for (i = 0; i < prime_limit; i++) 94 { 95 unsigned r; 96 /* FIXME: Reduce moduli + incr and store back, to allow for 97 division-free reductions. Alternatively, table primes[]'s 98 inverses (mod 2^16). */ 99 r = (moduli[i] + incr) % prime; 100 prime += primegap[i]; 101 102 if (r == 0) 103 goto next; 104 } 105 106 mpz_add_ui (p, p, difference); 107 difference = 0; 108 109 /* Miller-Rabin test */ 110 if (mpz_millerrabin (p, 25)) 111 goto done; 112 next:; 113 incr += 2; 114 } 115 mpz_add_ui (p, p, difference); 116 difference = 0; 117 } 118 done: 119 TMP_SFREE; 120} 121