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 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; } diff --git a/tests/mpz/t-nextprime.c b/tests/mpz/t-nextprime.c index 5607aea..4e8bfce 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,46 @@ 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 (); + } + + mpz_clear (y); + 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 (); } @@ -73,18 +128,43 @@ extern short diff4[]; extern short diff5[]; extern short diff6[]; -int -main (int argc, char **argv) +void +test_ref (gmp_randstate_ptr rands, int reps, + void (*func)(mpz_t, const mpz_t), + void(*ref_func)(mpz_t, const mpz_t)) { int i; - int reps = 20; - gmp_randstate_ptr rands; - mpz_t bs, x, nxtp, ref_nxtp; + mpz_t bs, x, test_p, ref_p; unsigned long size_range; - tests_start(); - rands = RANDS; + mpz_init (bs); + mpz_init (x); + mpz_init (test_p); + mpz_init (ref_p); + + 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) + abort (); + } + + mpz_clear (bs); + mpz_clear (x); + mpz_clear (test_p); + mpz_clear (ref_p); +} + +void +test_nextprime (gmp_randstate_ptr rands, int reps) +{ run ("2", 1000, "0x1ef7", diff1); run ("3", 1000 - 1, "0x1ef7", NULL); @@ -101,33 +181,70 @@ main (int argc, char **argv) run ("0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF80", 50, /* 2^128 - 128 */ "0x10000000000000000000000000000155B", diff6); - mpz_init (bs); - mpz_init (x); - mpz_init (nxtp); - mpz_init (ref_nxtp); + test_ref(rands, reps, mpz_nextprime, refmpz_nextprime); +} - TESTS_REPS (reps, argv, argc); +void +test_prevprime (gmp_randstate_ptr rands, int reps) +{ + int i, retval; + mpz_t x, prvp; - for (i = 0; i < reps; i++) - { - mpz_urandomb (bs, rands, 32); - size_range = mpz_get_ui (bs) % 8 + 2; /* 0..1024 bit operands */ + run_p ("2", 1000, "0x1ef7", diff1); - mpz_urandomb (bs, rands, size_range); - mpz_rrandomb (x, rands, mpz_get_ui (bs)); + run_p ("3", 1000 - 1, "0x1ef7", NULL); -/* gmp_printf ("%ld: %Zd\n", mpz_sizeinbase (x, 2), x); */ + run_p ("0x8a43866f5776ccd5b02186e90d28946aeb0ed914", 50, + "0x8a43866f5776ccd5b02186e90d28946aeb0eeec5", diff3); - mpz_nextprime (nxtp, x); - refmpz_nextprime (ref_nxtp, x); - if (mpz_cmp (nxtp, ref_nxtp) != 0) - abort (); + run_p ("0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF6C", 50, /* 2^148 - 148 */ + "0x100000000000000000000000000000000010ab", diff4); + + run_p ("0x1c2c26be55317530311facb648ea06b359b969715db83292ab8cf898d8b1b", 50, + "0x1c2c26be55317530311facb648ea06b359b969715db83292ab8cf898da957", diff5); + + run_p ("0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF80", 50, /* 2^128 - 128 */ + "0x10000000000000000000000000000155B", diff6); + + // Cast away int return from mpz_prevprime for test ref. + test_ref( + rands, reps, + (void (*)(mpz_t, const mpz_t)) mpz_prevprime, + refmpz_prevprime); + + // Test mpz_prevprime(x <= 2) returns 0, leaves rop unchanged. + mpz_init (x); + mpz_init (prvp); + mpz_set_ui (prvp, 123); + + 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 (); + } } - mpz_clear (bs); mpz_clear (x); - mpz_clear (nxtp); - mpz_clear (ref_nxtp); + mpz_clear (prvp); +} + +int +main (int argc, char **argv) +{ + gmp_randstate_ptr rands; + int reps = 20; + + tests_start(); + + rands = RANDS; + TESTS_REPS (reps, argv, argc); + + test_nextprime(rands, reps); + test_prevprime(rands, reps); tests_end (); return 0;