On Sun, Mar 15, 2020 at 4:38 PM Marco Bodrato <bodr...@mail.dm.unipi.it>
wrote:
> Ciao,
>
> Il 2020-03-15 00:07 Seth Troisi ha scritto:
> > New patch which is cleaned up and IMO ready to commit
>
> May I write a few more comments?
>
Always, my opinion about being ready is just that it's passed the bar of
being good enough not that it's perfect.
>
> > On Mon, Mar 9, 2020 at 5:15 PM Seth Troisi <brain...@gmail.com> wrote:
>
> >>>> I also dislike the boiler plate of the macros but I didn't
>
> Those macros should probably be moved to gmp-impl.h, to we avoid
> duplication.
> For this to be possible, we should avoid having different versions of
> the same macros in the different files...
>
I removed offset in the most recent version but I'm happy to add it back
and move the macros
into gmp-impl.h in a patch after this goes in.
>
> >> 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).
>
> Only the first "limb" of the sieve is hard-coded in primesieve.c, but it
> is possible to generate a larger "seed" at configure time, if we want.
>
I'm happy to look into this after this goes in.
>
> >>> Why not using a variable for INCR_LIMIT?
>
> I see now:
> + 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;
>
> Probably, we should change the name of the variable to
> odd_candidates_in_sieve, then you would probably have written:
> if (nbits <= 32)
> odd_numbers_in_sieve = 168;
> else if (nbits <= 64)
> odd_numbers_in_sieve = 775;
> else
> ...
>
Renamed to odds_in_composite_sieve
>
> Then I wonder if the cost of the /* Sieve next segment */ step is so
> high that we really have to try hard to avoid it. And if we are able to
> reduce the probability of its use to very unlikely, do we then really
> need to allocate the possibly huge next_mult array?
>
You're correct, this reduces the memory use 80%
>
> + mpz_root (tmp, tmp, 2);
>
> I'd use mpz_sqrt.
>
Done
>
> + // TODO: break the rare gap larger than this into gaps <= 255
>
> If you need to store larger gaps, you can also save in the array gap>>1,
> because they all are even. The prime limit will be able to go beyond
> 2^37...
>
Updated comment for how this could be done in the future.
>
> + /* Avoid having a primegap > 255, first occurence 436,273,009. */
> + ASSERT( 4000 < sieve_limit && sieve_limit < 436000000 );
>
> Why do you need to check (4000 < sieve_limit)?
>
It's a nice sanity check that sieve_limit was correctly calculated
>
> I'd expect a condition like
>
> + if (nbits <= numberof (primegap_small) * 2 + 1)
> + {
> + primegap = primegap_small;
> + prime_limit = nbits / 2;
> + }
>
> so that all primegap_small can be used. If you think this threshold is
> too high, then primegap_small can be reduced...
>
Done
Ĝis,
> m
>
diff -r bcca14c8a090 mpz/nextprime.c
--- a/mpz/nextprime.c Thu Mar 05 00:43:26 2020 +0100
+++ b/mpz/nextprime.c Sun Mar 15 18:23:48 2020 -0700
@@ -30,33 +30,109 @@
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,
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,
- 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,
- 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,
- 6,14,4,6,6,8,6,12
+ 12,2,18,6,10
};
-#define NUMBER_OF_PRIMES 167
+#define NUMBER_OF_PRIMES 100
+
+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_sqrt (tmp, tmp);
+ mpz_tdiv_q_ui(tmp, tmp, 124);
+
+ // TODO: Storing gap/2 would allow this to go above 436M.
+ 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( 1000 < 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;
+ 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 odds_in_composite_sieve;
+ unsigned long difference;
+ TMP_DECL;
/* First handle tiny numbers */
if (mpz_cmp_ui (n, 2) < 0)
@@ -70,59 +146,108 @@
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 (2 * nbits <= NUMBER_OF_PRIMES)
+ {
+ 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;
- TMP_SMARK;
+ /* Both primesieve and prime_limit ignore the first two primes. */
+ ASSERT(i == prime_limit);
+ }
+
+ if (nbits <= 32)
+ odds_in_composite_sieve = 336 / 2 + 1;
+ else if (nbits <= 64)
+ odds_in_composite_sieve = 1550 / 2 + 1;
+ else
+ /* Corresponds to a merit 14 prime_gap, which is rare. */
+ odds_in_composite_sieve = 5 * nbits;
+
+ composite = TMP_SALLOC_TYPE (odds_in_composite_sieve, char);
+
+ /* composite[2*i] stores if p+2*i is a known composite */
+ memset (composite, 0, odds_in_composite_sieve);
- /* Compute residues modulo small odd primes */
- moduli = TMP_SALLOC_TYPE (prime_limit, unsigned short);
+ prime = 3;
+ for (i = 0; i < prime_limit; i++)
+ {
+ m = mpz_cdiv_ui(p, prime);
+ /* Only care about odd multiplies of prime. */
+ if (m & 1)
+ m += prime;
+ m >>= 1;
+ /* mark off any composites in sieve */
+ for (; m < odds_in_composite_sieve; m += prime)
+ composite[m] = 1;
+ 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];
- }
-
-#define INCR_LIMIT 0x10000 /* deep science */
+ for (incr = 0; incr < odds_in_composite_sieve; difference += 2, incr += 1)
+ {
+ if (composite[incr])
+ continue;
- 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];
+ mpz_add_ui (p, p, difference);
+ difference = 0;
- if (r == 0)
- goto next;
- }
-
- mpz_add_ui (p, p, difference);
- difference = 0;
+ /* Miller-Rabin test */
+ if (mpz_millerrabin (p, 25))
+ goto done;
+ }
- /* Miller-Rabin test */
- if (mpz_millerrabin (p, 25))
- goto done;
- next:;
- incr += 2;
- }
+ /* Sieve next segment, very rare */
mpz_add_ui (p, p, difference);
difference = 0;
+
+ prime = 3;
+ memset (composite, 0, odds_in_composite_sieve);
+ 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 < odds_in_composite_sieve; m += prime)
+ composite[m] = 1;
+ 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
