I wrote up an initial implementation, tests, and documentation https://github.com/sethtroisi/libgmp/pull/10/files
On Tue, Sep 3, 2019 at 3:19 AM Pedro Gimeno <gmpde...@formauri.es> wrote: > Seth Troisi wrote on 2019-09-03 08:59: > > > 3. add mpz_nextprime_seq(a, b) => next prime a + k*b with k >= 1 > > Personally I would prefer it named mpz_nextprime_step or > mpz_nextprime_stride or mpz_nextprime_strided. Or even > mpz_nexprime_with_step, if that doesn't violate some naming convention I'm > not aware of. It's not all that clear to me what sequence does 'seq' refer > to in the name. >
From 353f07c97a7fa142c32119c87ffba0fe925e70bb Mon Sep 17 00:00:00 2001 From: Seth Troisi <sethtro...@google.com> Date: Wed, 4 Sep 2019 23:46:10 -0700 Subject: [PATCH 1/3] Add mpz_prevprime --- gmp-h.in | 3 ++ mpz/nextprime.c | 129 +++++++++++++++++++++++++++++++++--------------- 2 files changed, 93 insertions(+), 39 deletions(-) diff --git a/gmp-h.in b/gmp-h.in index f448b4e..a1208d8 100644 --- a/gmp-h.in +++ b/gmp-h.in @@ -947,6 +947,9 @@ __GMP_DECLSPEC void mpz_neg (mpz_ptr, mpz_srcptr); #define mpz_nextprime __gmpz_nextprime __GMP_DECLSPEC void mpz_nextprime (mpz_ptr, mpz_srcptr); +#define mpz_prevprime __gmpz_prevprime +__GMP_DECLSPEC int mpz_prevprime (mpz_ptr, mpz_srcptr); + #define mpz_out_raw __gmpz_out_raw #ifdef _GMP_H_HAVE_FILE __GMP_DECLSPEC size_t mpz_out_raw (FILE *, mpz_srcptr); diff --git a/mpz/nextprime.c b/mpz/nextprime.c index 4c5ca57..90f03a4 100644 --- a/mpz/nextprime.c +++ b/mpz/nextprime.c @@ -33,6 +33,49 @@ see https://www.gnu.org/licenses/. */ #include "gmp-impl.h" #include "longlong.h" +#define INCR_LIMIT 0x10000 /* deep science */ + +/* Loop to find next/prev prime */ +#define next_prime_helper(op, function) \ + do { \ + /* compute residues modulo small odd primes */ \ + moduli = TMP_SALLOC_TYPE (prime_limit, unsigned short); \ + \ + for (;;) \ + { \ + prime = 3; \ + for (i = 0; i < prime_limit; i++) \ + { \ + moduli[i] = mpz_tdiv_ui (p, prime); \ + prime += primegap[i]; \ + } \ + for (difference = incr = 0; incr < INCR_LIMIT; difference += 2) \ + { \ + /* First check residues */ \ + prime = 3; \ + for (i = 0; i < prime_limit; i++) \ + { \ + signed r; \ + signed t = moduli[i] op incr; \ + r = t % ((int) prime); \ + prime += primegap[i]; \ + if (r == 0) \ + goto next; \ + } \ + function (p, p, difference); \ + difference = 0; \ + \ + /* Miller-Rabin test */ \ + if (mpz_millerrabin (p, 25)) \ + goto done; \ + next:; \ + incr += 2; \ + } \ + function (p, p, difference); \ + difference = 0; \ + } \ + } while (0) + static const unsigned char primegap[] = { 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, @@ -79,50 +122,58 @@ mpz_nextprime (mpz_ptr p, mpz_srcptr n) TMP_SMARK; - /* Compute residues modulo small odd primes */ - moduli = TMP_SALLOC_TYPE (prime_limit * sizeof moduli[0], unsigned short); + /* work forwards looking for number with moduli == 0 and prime */ + next_prime_helper (+, mpz_add_ui); + + done: + TMP_SFREE; +} + +int +mpz_prevprime (mpz_ptr p, mpz_srcptr n) +{ + unsigned short *moduli; + unsigned long difference; + int i; + unsigned prime_limit; + unsigned long prime; + mp_size_t pn; + mp_bitcnt_t nbits; + unsigned incr; + TMP_SDECL; + + /* handle numbers with no previous prime. */ + if (mpz_cmp_ui (n, 2) <= 0) + return 0; - for (;;) + if (mpz_cmp_ui (n, 3) <= 0) { - /* FIXME: Compute lazily? */ - prime = 3; - for (i = 0; i < prime_limit; i++) - { - moduli[i] = mpz_tdiv_ui (p, prime); - prime += primegap[i]; - } + /* previous prime for 3 == 2 */ + mpz_set_ui (p, 2); + return 2; + } -#define INCR_LIMIT 0x10000 /* deep science */ + /* First odd less than n */ + mpz_sub_ui (p, n, 2); + mpz_setbit (p, 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]; - - 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; + if (mpz_cmp_ui (p, 7) <= 0) + { + return 2; } + + pn = SIZ(p); + MPN_SIZEINBASE_2EXP(nbits, PTR(p), pn, 1); + if (nbits / 2 >= NUMBER_OF_PRIMES) + prime_limit = NUMBER_OF_PRIMES - 1; + else + prime_limit = nbits / 2; + + TMP_SMARK; + + /* work backwards looking for number with moduli == 0 and prime */ + next_prime_helper (-, mpz_sub_ui); + done: TMP_SFREE; } From 15bccec04a7f4f43021f59f247ee2fc2a9e43d2a Mon Sep 17 00:00:00 2001 From: Seth Troisi <sethtro...@google.com> Date: Wed, 4 Sep 2019 23:48:19 -0700 Subject: [PATCH 2/3] Added tests --- tests/mpz/t-nextprime.c | 245 +++++++++++++++++++++++++++++++++++----- 1 file changed, 215 insertions(+), 30 deletions(-) diff --git a/tests/mpz/t-nextprime.c b/tests/mpz/t-nextprime.c index 5607aea..0f82be9 100644 --- a/tests/mpz/t-nextprime.c +++ b/tests/mpz/t-nextprime.c @@ -32,6 +32,23 @@ refmpz_nextprime (mpz_ptr p, mpz_srcptr t) mpz_add_ui (p, p, 1L); } +void +refmpz_prevprime (mpz_ptr p, mpz_srcptr t) +{ + if (mpz_cmp_ui(t, 2) <= 0) + return; + + if (mpz_cmp_ui(t, 3) <= 0) + { + mpz_set_ui (p, 2); + return; + } + + mpz_sub_ui (p, t, 1L); + while (! mpz_probab_prime_p (p, 10)) + mpz_sub_ui (p, p, 1L); +} + void run (const char *start, int reps, const char *end, short diffs[]) { @@ -58,8 +75,8 @@ run (const char *start, int reps, const char *end, short diffs[]) if (mpz_cmp (x, y) != 0) { - gmp_printf ("got %Zx\n", x); - gmp_printf ("want %Zx\n", y); + gmp_printf ("got %Zd\n", x); + gmp_printf ("want %Zd\n", y); abort (); } @@ -67,24 +84,87 @@ run (const char *start, int reps, const char *end, short diffs[]) mpz_clear (x); } +void +run_p (const char *start, int reps, const char *end, short diffs[]) +{ + mpz_t x, y; + int i; + + mpz_init_set_str (x, end, 0); + mpz_init (y); + + // Last rep doesn't share same data with nextprime + for (i = 0; i < reps - 1; i++) + { + mpz_prevprime (y, x); + mpz_sub (x, x, y); + if (diffs != NULL && + (! mpz_fits_sshort_p (x) || diffs[reps - i - 1] != (short) mpz_get_ui (x))) + { + gmp_printf ("diff list discrepancy\n"); + abort (); + } + mpz_swap (x, y); + } + + // starts aren't always prime, so check that result is less than or equal + mpz_prevprime(x, x); + + mpz_set_str(y, start, 0); + if (mpz_cmp (x, y) > 0) + { + gmp_printf ("got %Zd\n", x); + gmp_printf ("want %Zd\n", y); + abort (); + } + + mpz_clear (y); + mpz_clear (x); +} + +void +test_ref (gmp_randstate_ptr rands, int reps, void (*func)(mpz_t, mpz_t), void(*ref_func)(mpz_t, mpz_t)) +{ + int i; + mpz_t bs, x, test_p, ref_p; + unsigned long size_range; + + mpz_inits (bs, x, test_p, ref_p, NULL); + + mpz_set_ui(test_p, 0); + mpz_set_ui(ref_p, 0); + + for (i = 0; i < reps; i++) + { + mpz_urandomb (bs, rands, 32); + size_range = mpz_get_ui (bs) % 8 + 2; /* 0..1024 bit operands */ + + mpz_urandomb (bs, rands, size_range); + mpz_rrandomb (x, rands, mpz_get_ui (bs)); + + func (test_p, x); + ref_func (ref_p, x); + if (mpz_cmp (test_p, ref_p) != 0) + { + gmp_printf ("op %Zd\n", x); + gmp_printf ("got %Zd\n", test_p); + gmp_printf ("want %Zd\n", ref_p); + abort (); + } + } + + mpz_clears (bs, x, test_p, ref_p, NULL); +} + extern short diff1[]; extern short diff3[]; extern short diff4[]; extern short diff5[]; extern short diff6[]; -int -main (int argc, char **argv) +void +test_nextprime (gmp_randstate_ptr rands, int reps) { - int i; - int reps = 20; - gmp_randstate_ptr rands; - mpz_t bs, x, nxtp, ref_nxtp; - unsigned long size_range; - - tests_start(); - rands = RANDS; - run ("2", 1000, "0x1ef7", diff1); run ("3", 1000 - 1, "0x1ef7", NULL); @@ -101,33 +181,138 @@ main (int argc, char **argv) run ("0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF80", 50, /* 2^128 - 128 */ "0x10000000000000000000000000000155B", diff6); - mpz_init (bs); + test_ref(rands, reps, mpz_nextprime, refmpz_nextprime); +} + +void +test_prevprime (gmp_randstate_ptr rands, int reps) +{ + int i, retval; + mpz_t x, prvp; + + run_p ("2", 1000, "0x1ef7", diff1); + + run_p ("3", 1000 - 1, "0x1ef7", NULL); + + run_p ("0x8a43866f5776ccd5b02186e90d28946aeb0ed914", 50, + "0x8a43866f5776ccd5b02186e90d28946aeb0eeec5", diff3); + + run_p ("0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF6C", 50, /* 2^148 - 148 */ + "0x100000000000000000000000000000000010ab", diff4); + + run_p ("0x1c2c26be55317530311facb648ea06b359b969715db83292ab8cf898d8b1b", 50, + "0x1c2c26be55317530311facb648ea06b359b969715db83292ab8cf898da957", diff5); + + run_p ("0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF80", 50, /* 2^128 - 128 */ + "0x10000000000000000000000000000155B", diff6); + + test_ref(rands, reps, mpz_prevprime, refmpz_prevprime); + + // Test mpz_prevprime(x <= 2) returns 0, leaves rop unchanged. mpz_init (x); - mpz_init (nxtp); - mpz_init (ref_nxtp); + mpz_init (prvp); + mpz_set_ui (prvp, 123); - TESTS_REPS (reps, argv, argc); + for (i = -10; i <= 2; i++) + { + mpz_set_si(x, i); + retval = mpz_prevprime (prvp, x); + if ( retval != 0 || mpz_get_si (prvp) != 123 ) + { + gmp_printf ("mpz_prevprime(%Zd) return (%d) rop (%Zd)\n", x, prvp ); + abort (); + } + } - for (i = 0; i < reps; i++) + mpz_clear (x); + mpz_clear (prvp); +} + +void +test_largegap (mpz_t low, const gap) +{ + mpz_t t, nxt, prv; + + mpz_inits (t, nxt, prv, NULL); + + mpz_nextprime(nxt, low); + mpz_sub(t, nxt, low); + + if (mpz_cmp_ui(t, gap) != 0) { - mpz_urandomb (bs, rands, 32); - size_range = mpz_get_ui (bs) % 8 + 2; /* 0..1024 bit operands */ + gmp_printf ("prime gap %Zd != %d\n", t, gap); + abort (); + } - mpz_urandomb (bs, rands, size_range); - mpz_rrandomb (x, rands, mpz_get_ui (bs)); + mpz_prevprime(prv, nxt); + + if (mpz_cmp(prv, low)) + { + gmp_printf ("mpz_prevprime(mpz_nextprime(x)) != x\n"); + abort (); + } + + mpz_clears (t, nxt, prv, NULL); +} + +void +test_largegaps () +{ + mpz_t x, t; + + mpz_init (x); + mpz_init (t); -/* gmp_printf ("%ld: %Zd\n", mpz_sizeinbase (x, 2), x); */ + /* + // This takes ~90 seconds, it test the deep science magic constant in + // nextprime.c but takes too long to be always enabled. - mpz_nextprime (nxtp, x); - refmpz_nextprime (ref_nxtp, x); - if (mpz_cmp (nxtp, ref_nxtp) != 0) - abort (); + // Gap 66520 from P816 = 1931 * 1933# / 7230 - 30244 + mpz_set_ui (x, 1931); + mpz_set_ui (t, 2); + while (mpz_cmp_ui(t, 1933) <= 0) + { + mpz_mul (x, x, t); + mpz_nextprime (t, t); } + mpz_divexact_ui (x, x, 7230); + mpz_sub_ui (x, x, 30244); + + test_largegap(x, 66520); + */ + + // Gap 33008 from P454 = 55261931 * 1063#/210 - 13116 + mpz_set_ui (x, 55261931); + mpz_set_ui (t, 2); + while (mpz_cmp_ui(t, 1063) <= 0) + { + mpz_mul (x, x, t); + mpz_nextprime (t, t); + } + mpz_divexact_ui (x, x, 210); + mpz_sub_ui (x, x, 13116); + + test_largegap(x, 33008); + - mpz_clear (bs); mpz_clear (x); - mpz_clear (nxtp); - mpz_clear (ref_nxtp); + mpz_clear (t); +} +int +main (int argc, char **argv) +{ + gmp_randstate_ptr rands; + int reps = 500; + + tests_start(); + + rands = RANDS; + TESTS_REPS (reps, argv, argc); + + test_nextprime(rands, reps); + test_prevprime(rands, reps); + + test_largegaps(); tests_end (); return 0; From 7f563a4494038d6efb7312c16550423a3097f5e6 Mon Sep 17 00:00:00 2001 From: Seth Troisi <sethtro...@google.com> Date: Wed, 4 Sep 2019 23:48:33 -0700 Subject: [PATCH 3/3] Added documentation --- doc/gmp.texi | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/doc/gmp.texi b/doc/gmp.texi index 60c2634..da412a1 100644 --- a/doc/gmp.texi +++ b/doc/gmp.texi @@ -3559,8 +3559,19 @@ and 50. @deftypefun void mpz_nextprime (mpz_t @var{rop}, const mpz_t @var{op}) @cindex Next prime function Set @var{rop} to the next prime greater than @var{op}. +@end deftypefun + +@deftypefun int mpz_prevprime (mpz_t @var{rop}, const mpz_t @var{op}) +@cindex Previous prime function +Set @var{rop} to the greatest prime less than @var{op}. + +If previous prime doesn't exist (e.g. @var{op} < 3), rop is unchanged and +0 is returned. + +Return 1 if @var{rop} is a probably prime, and 2 if @var{rop} is definitely +prime. -This function uses a probabilistic algorithm to identify primes. For +These functions use a probabilistic algorithm to identify primes. For practical purposes it's adequate, the chance of a composite passing will be extremely small. @end deftypefun
_______________________________________________ gmp-devel mailing list gmp-devel@gmplib.org https://gmplib.org/mailman/listinfo/gmp-devel
