Author: cperciva
Date: Sun Jun 4 02:36:37 2017
New Revision: 319561
URL: https://svnweb.freebsd.org/changeset/base/319561
Log:
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
Modified:
head/usr.bin/primes/primes.c
head/usr.bin/primes/primes.h
head/usr.bin/primes/spsp.c
Modified: head/usr.bin/primes/primes.c
==============================================================================
--- head/usr.bin/primes/primes.c Sun Jun 4 02:21:38 2017 (r319560)
+++ head/usr.bin/primes/primes.c Sun Jun 4 02:36:37 2017 (r319561)
@@ -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. */
Modified: head/usr.bin/primes/primes.h
==============================================================================
--- head/usr.bin/primes/primes.h Sun Jun 4 02:21:38 2017 (r319560)
+++ head/usr.bin/primes/primes.h Sun Jun 4 02:36:37 2017 (r319561)
@@ -72,6 +72,3 @@ extern const size_t pattern_size; /* length of pattern
/* Test for primality using strong pseudoprime tests. */
int isprime(ubig);
-
-/* Maximum value which the SPSP code can handle. */
-#define SPSPMAX 3825123056546413050ULL
Modified: head/usr.bin/primes/spsp.c
==============================================================================
--- head/usr.bin/primes/spsp.c Sun Jun 4 02:21:38 2017 (r319560)
+++ head/usr.bin/primes/spsp.c Sun Jun 4 02:36:37 2017 (r319561)
@@ -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);
}