New patch which is cleaned up and IMO ready to commit On Mon, Mar 9, 2020 at 5:15 PM Seth Troisi <brain...@gmail.com> wrote:
> > > On Mon, Mar 9, 2020 at 5:03 AM Marco Bodrato <bodr...@mail.dm.unipi.it> > wrote: > >> Ciao, >> >> Il 2020-03-09 11:59 Seth Troisi ha scritto: >> > On Mon, Mar 9, 2020 at 2:05 AM Marco Bodrato >> > <bodr...@mail.dm.unipi.it> wrote: >> > >> >> Ciao, >> >> >> >> Il 2020-03-09 02:56 Seth Troisi ha scritto: >> >> > It's highly possible I misunderstand how these macros are supposed to >> > be used. >> > I also dislike the boiler plate of the macros but I didn't like >> >> When I say "that code is messy", I mean my code. >> And you agree :-) those macros are just easy to misunderstand :-/ >> >> > If you have suggestions or code pointers for a better pattern I'd >> > appreciate that. >> >> >> Did you try to use the gmp_primesieve function directly? >> > >> > I'm not sure what you're referencing here, I tried with >> >> I'm proposing to get rid of primegap, and use something like the >> following. >> This is just a copy-paste from your code, not a really working sequence! >> >> + __sieve = TMP_ALLOC_LIMBS (primesieve_size (sieve_limit)); >> + prime_limit_pre = gmp_primesieve(__sieve, sieve_limit); >> >> + next_mult = TMP_ALLOC_TYPE (prime_limit, unsigned int); >> [...] >> + /* Invert p for modulo */ >> + mpz_neg(p, p); >> >> [..handle 3 outside the loop] >> >> + LOOP_ON_SIEVE_BEGIN (prime, n_to_bit (5), n_to_bit (sieve_limit), >> 0, __sieve); >> + m = mpz_fdiv_ui(p, prime); >> + /* Only care about odd multiplies (even indexes) */ >> + if (m & 1) >> + m += prime; >> + >> + /* mark off any composites in sieve */ >> + for (; m < INCR_LIMIT; m += 2*prime) >> + sieve[m] = 1; >> + next_mult[i] = m; >> + LOOP_ON_SIEVE_END; >> + mpz_neg(p, p); >> >> [...] >> >> + /* Sieve next segment */ >> [..handle 3] >> + memset(sieve, 0, INCR_LIMIT); >> + LOOP_ON_SIEVE_BEGIN (prime, n_to_bit (5), n_to_bit (sieve_limit), >> 0, __sieve); >> + m = next_mult[i] - INCR_LIMIT; >> + for (; m < INCR_LIMIT; m += prime * 2) >> + sieve[m] = 1; >> + next_mult[i] = m; >> + LOOP_ON_SIEVE_END; >> > When nbits is small (IMHO the most common case) we'd like to be able to > fallback to a constant array which I don't think exist for primesieve > (because gmp_limb_bits isn't constant). > So I lowered the threshold to 150M (50MB) in the future I have a trick > that will reduce memory usage and we can try to increase this limit. > > >> > Is there a way to unalloc a TMP_ALLOC_LIMBS before TMP_FREE? >> >> No, if you want to unalloc, you must use __GMP_ALLOCATE_FUNC_TYPE, >> __GMP_REALLOCATE_FUNC_TYPE, and __GMP_FREE_FUNC_TYPE. >> >> > Sieve next segment is rare, the average gap is ln(n), if INCR_LIMIT >> > was changed to be a variable it could be set >> > so /*sieve next segment*/ happened on average less than once, and then >> >> Why not using a variable for INCR_LIMIT? >> > It will be another change with more moving parts. I'm hoping to do it > after the large part of this change goes in. > >> >> >> In the array char *sieve, you use only the even entries. Maybe you >> >> can reduce its size by half, and remove some "*2" around the code? >> > >> > it's only 4000 entries so I wasn't bothering at this time, but it >> >> I agree. >> > Fixed now. > >> >> PS: did you consider mpz_Cdiv_ui, instead of _neg,_Fdiv_ui,_neg ? >> > Nope, this is much cleaner. Thanks. > >> >> Ĝis, >> m >> >
diff -r bcca14c8a090 mpz/nextprime.c --- a/mpz/nextprime.c Thu Mar 05 00:43:26 2020 +0100 +++ b/mpz/nextprime.c Sat Mar 14 15:56:07 2020 -0700 @@ -30,10 +30,52 @@ GNU Lesser General Public License along with the GNU MP Library. If not, see https://www.gnu.org/licenses/. */ +#include <string.h> + #include "gmp-impl.h" #include "longlong.h" -static const unsigned char primegap[] = +/*********************************************************/ +/* Section sieve: sieving functions and tools for primes */ +/*********************************************************/ + +static mp_limb_t +id_to_n (mp_limb_t id) { return id*3+1+(id&1); } + +static mp_limb_t +n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; } + +static mp_size_t +primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; } + + +/************************************/ +/* Section macros: macros for sieve */ +/************************************/ + +#define LOOP_ON_SIEVE_BEGIN(prime,start,end,sieve) \ + do { \ + mp_limb_t __mask, __index, __max_i, __i; \ + __i = (start); \ + __index = __i / GMP_LIMB_BITS; \ + __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS); \ + __max_i = (end); \ + do { \ + ++__i; \ + if (((sieve)[__index] & __mask) == 0) \ + { \ + mp_limb_t prime = id_to_n(__i) \ + +#define LOOP_ON_SIEVE_END \ + } \ + __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1); \ + __index += __mask & 1; \ + } while (__i <= __max_i); \ + } while (0) + + + +static const unsigned char primegap_small[] = { 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, 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, @@ -45,18 +87,55 @@ #define NUMBER_OF_PRIMES 167 +long calculate_sievelimit(mp_bitcnt_t nbits) { + long sieve_limit; + + // Estimate a good sieve bound. Based on derivative of + // Merten's 3rd theorem * avg gap * cost of mod + // vs + // Cost of PRP test O(N^2.55) + if (nbits < 12800) + { + mpz_t tmp; + // sieve_limit ~= nbits ^ (5/2) / 124 + mpz_init (tmp); + mpz_ui_pow_ui (tmp, nbits, 5); + mpz_root (tmp, tmp, 2); + mpz_tdiv_q_ui(tmp, tmp, 124); + + // TODO: break the rare gap larger than this into gaps <= 255 + sieve_limit = mpz_get_ui(tmp); + mpz_clear (tmp); + } + else + { + /* Larger threshold is faster but takes ~O(5*n/ln(n)) memory. + * For 33,000 bits limitting to 150M is ~12% slower than using the + * optimal 1.5B sieve_limit. + */ + sieve_limit = 150000000; + } + + /* Avoid having a primegap > 255, first occurence 436,273,009. */ + ASSERT( 4000 < sieve_limit && sieve_limit < 436000000 ); + return sieve_limit; +} + void mpz_nextprime (mpz_ptr p, mpz_srcptr n) { - unsigned short *moduli; - unsigned long difference; - int i; + unsigned int *next_mult; + char *composite; + const unsigned char *primegap; unsigned prime_limit; - unsigned long prime; + unsigned prime; mp_size_t pn; mp_bitcnt_t nbits; + int i, m; unsigned incr; - TMP_SDECL; + unsigned incr_limit; + unsigned long difference; + TMP_DECL; /* First handle tiny numbers */ if (mpz_cmp_ui (n, 2) < 0) @@ -70,59 +149,105 @@ if (mpz_cmp_ui (p, 7) <= 0) return; + TMP_MARK; pn = SIZ(p); MPN_SIZEINBASE_2EXP(nbits, PTR(p), pn, 1); - if (nbits / 2 >= NUMBER_OF_PRIMES) - prime_limit = NUMBER_OF_PRIMES - 1; + if (nbits <= 199) + { + primegap = primegap_small; + prime_limit = nbits / 2; + } else - prime_limit = nbits / 2; + { + long sieve_limit; + mp_limb_t *sieve; + unsigned char *primegap_tmp; + int last_prime; + + /* sieve numbers up to sieve_limit and save prime count */ + sieve_limit = calculate_sievelimit(nbits); + sieve = TMP_ALLOC_LIMBS (primesieve_size (sieve_limit)); + prime_limit = gmp_primesieve(sieve, sieve_limit); + + /* Needed to avoid assignment of read-only location */ + primegap_tmp = TMP_ALLOC_TYPE (prime_limit, unsigned char); + primegap = primegap_tmp; + + i = 0; + last_prime = 3; + LOOP_ON_SIEVE_BEGIN (prime, n_to_bit (5), n_to_bit (sieve_limit), sieve); + primegap_tmp[i++] = prime - last_prime; + last_prime = prime; + LOOP_ON_SIEVE_END; + + /* Both primesieve and prime_limit ignore the first two primes. */ + ASSERT(i == prime_limit); + } - TMP_SMARK; + /* Compute next multiple for small odd primes */ + next_mult = TMP_ALLOC_TYPE (prime_limit, unsigned int); + + if (nbits <= 32) + incr_limit = 336; + else if (nbits <= 64) + incr_limit = 1550; + else + /* Corresponds to a merit 10 prime_gap, which is rare. */ + incr_limit = 7 * nbits / 2; + + composite = TMP_SALLOC_TYPE (incr_limit, char); + + /* composite[2*i] remembers if p+2*i is a known composite */ + memset (composite, 0, incr_limit); - /* Compute residues modulo small odd primes */ - moduli = TMP_SALLOC_TYPE (prime_limit, unsigned short); + /* Invert p for modulo */ + prime = 3; + for (i = 0; i < prime_limit; i++) + { + m = mpz_cdiv_ui(p, prime); + /* Only care about odd multiplies */ + if (m & 1) + m += prime; + m >>= 1; + /* mark off any composites in sieve */ + for (; m < incr_limit; m += prime) + composite[m] = 1; + next_mult[i] = m; + prime += primegap[i]; + } + + difference = 0; for (;;) { - /* FIXME: Compute lazily? */ - prime = 3; - for (i = 0; i < prime_limit; i++) - { - moduli[i] = mpz_tdiv_ui (p, prime); - prime += primegap[i]; - } + for (incr = 0; incr < incr_limit; difference += 2, incr += 1) + { + if (composite[incr]) + continue; -#define INCR_LIMIT 0x10000 /* deep science */ + mpz_add_ui (p, p, difference); + difference = 0; - for (difference = incr = 0; incr < INCR_LIMIT; difference += 2) - { - /* First check residues */ - prime = 3; - for (i = 0; i < prime_limit; i++) - { - unsigned r; - /* FIXME: Reduce moduli + incr and store back, to allow for - division-free reductions. Alternatively, table primes[]'s - inverses (mod 2^16). */ - r = (moduli[i] + incr) % prime; - prime += primegap[i]; + /* Miller-Rabin test */ + if (mpz_millerrabin (p, 25)) + goto done; + } - if (r == 0) - goto next; - } - - mpz_add_ui (p, p, difference); - difference = 0; - - /* Miller-Rabin test */ - if (mpz_millerrabin (p, 25)) - goto done; - next:; - incr += 2; - } - mpz_add_ui (p, p, difference); - difference = 0; + /* Sieve next segment */ + prime = 3; + memset (composite, 0, incr_limit); + for (i = 0; i < prime_limit; i++) + { + m = next_mult[i] - incr_limit; + for (; m < incr_limit; m += prime) + composite[m] = 1; + next_mult[i] = m; + prime += primegap[i]; + } } done: - TMP_SFREE; + TMP_FREE; } + +#undef LOOP_ON_SIEVE_END +#undef LOOP_ON_SIEVE_BEGIN
_______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel
