Using results from

J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime
    bases, Math. Comp. 86(304):985-1003, 2017.
teach primes(6) to enumerate primes up to 2^64 - 1.  Until Sorenson
and Webster's paper, we did not know how many strong speudoprime tests
were required when testing alleged primes between 3825123056546413051
and 2^64 - 1.

Reported by:	Luiz Henrique de Figueiredo
Relnotes:	primes(6) now enumerates primes beyond 3825123056546413050,
		up to a new limit of 2^64 - 1.
MFC After:	1 week
This commit is contained in:
Colin Percival 2017-06-04 02:36:37 +00:00
parent 86f2d54ff0
commit ade8bcee50
3 changed files with 30 additions and 15 deletions

View File

@ -120,7 +120,7 @@ main(int argc, char *argv[])
argv += optind;
start = 0;
stop = SPSPMAX;
stop = (uint64_t)(-1);
/*
* Convert low and high args. Strtoumax(3) sets errno to
@ -147,8 +147,6 @@ main(int argc, char *argv[])
err(1, "%s", argv[1]);
if (*p != '\0')
errx(1, "%s: illegal numeric format.", argv[1]);
if (stop > SPSPMAX)
errx(1, "%s: stop value too large.", argv[1]);
break;
case 1:
/* Start on the command line. */

View File

@ -72,6 +72,3 @@ extern const size_t pattern_size; /* length of pattern array */
/* Test for primality using strong pseudoprime tests. */
int isprime(ubig);
/* Maximum value which the SPSP code can handle. */
#define SPSPMAX 3825123056546413050ULL

View File

@ -32,23 +32,32 @@ __FBSDID("$FreeBSD$");
#include "primes.h"
/* Return a * b % n, where 0 <= a, b < 2^63, 0 < n < 2^63. */
/* Return a * b % n, where 0 < n. */
static uint64_t
mulmod(uint64_t a, uint64_t b, uint64_t n)
{
uint64_t x = 0;
uint64_t an = a % n;
while (b != 0) {
if (b & 1)
x = (x + a) % n;
a = (a + a) % n;
if (b & 1) {
x += an;
if ((x < an) || (x >= n))
x -= n;
}
if (an + an < an)
an = an + an - n;
else if (an + an >= n)
an = an + an - n;
else
an = an + an;
b >>= 1;
}
return (x);
}
/* Return a^r % n, where 0 <= a < 2^63, 0 < n < 2^63. */
/* Return a^r % n, where 0 < n. */
static uint64_t
powmod(uint64_t a, uint64_t r, uint64_t n)
{
@ -173,9 +182,20 @@ isprime(ubig _n)
if (n < 3825123056546413051)
return (1);
/* We can't handle values larger than this. */
assert(n <= SPSPMAX);
/*
* Value from:
* J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime
* bases, Math. Comp. 86(304):985-1003, 2017.
*/
/* UNREACHABLE */
return (0);
/* No SPSPs to bases 2..37 less than 318665857834031151167461. */
if (!spsp(n, 29))
return (0);
if (!spsp(n, 31))
return (0);
if (!spsp(n, 37))
return (0);
/* All 64-bit values are less than 318665857834031151167461. */
return (1);
}