summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorcperciva <cperciva@FreeBSD.org>2017-06-22 05:26:08 +0000
committercperciva <cperciva@FreeBSD.org>2017-06-22 05:26:08 +0000
commit17d16f2e819d8dfb24936b38c0f86c0a399340be (patch)
tree8866bb504ec4fc4f99e80cd102b1d5d80cd88e98
parentdc58c32e0b958b3dae4307ccd29d677a51ce3cb5 (diff)
downloadFreeBSD-src-17d16f2e819d8dfb24936b38c0f86c0a399340be.zip
FreeBSD-src-17d16f2e819d8dfb24936b38c0f86c0a399340be.tar.gz
MFC r31956[12]: Teach primes(6) to enumerate primes up to 2^64 - 1.
Approved by: re (delphij) Relnotes: primes(6) now enumerates primes beyond 3825123056546413050, up to a new limit of 2^64 - 1. > Description of fields to fill in above: 76 columns --| > PR: If and which Problem Report is related. > Submitted by: If someone else sent in the change. > Reported by: If someone else reported the issue. > Reviewed by: If someone else reviewed your modification. > Approved by: If you needed approval for this commit. > Obtained from: If the change is from a third party. > MFC after: N [day[s]|week[s]|month[s]]. Request a reminder email. > MFH: Ports tree branch name. Request approval for merge. > Relnotes: Set to 'yes' for mention in release notes. > Security: Vulnerability reference (one per line) or description. > Sponsored by: If the change was sponsored by an organization. > Differential Revision: https://reviews.freebsd.org/D### (*full* phabric URL needed). > Empty fields above will be automatically removed. _M . M usr.bin/factor/factor.6 M usr.bin/primes/primes.c M usr.bin/primes/primes.h M usr.bin/primes/spsp.c
-rw-r--r--usr.bin/factor/factor.66
-rw-r--r--usr.bin/primes/primes.c4
-rw-r--r--usr.bin/primes/primes.h3
-rw-r--r--usr.bin/primes/spsp.c38
4 files changed, 30 insertions, 21 deletions
diff --git a/usr.bin/factor/factor.6 b/usr.bin/factor/factor.6
index ba82f14..6d01a17 100644
--- a/usr.bin/factor/factor.6
+++ b/usr.bin/factor/factor.6
@@ -119,9 +119,3 @@ cannot handle the
factor list,
.Nm primes
will not get you a world record.
-.Pp
-.Nm primes
-is unable to list primes between 3825123056546413050 and 18446744073709551615
-since it relies on strong pseudoprime tests after sieving, and nobody has
-proven how many strong pseudoprime tests are required to prove primality for
-integers larger than 3825123056546413050.
diff --git a/usr.bin/primes/primes.c b/usr.bin/primes/primes.c
index a1c95c2..3a6e85f 100644
--- a/usr.bin/primes/primes.c
+++ b/usr.bin/primes/primes.c
@@ -112,7 +112,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
@@ -139,8 +139,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. */
diff --git a/usr.bin/primes/primes.h b/usr.bin/primes/primes.h
index 3a18fc7..ae0bbcb 100644
--- a/usr.bin/primes/primes.h
+++ b/usr.bin/primes/primes.h
@@ -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
diff --git a/usr.bin/primes/spsp.c b/usr.bin/primes/spsp.c
index f61acd6..117333b 100644
--- a/usr.bin/primes/spsp.c
+++ b/usr.bin/primes/spsp.c
@@ -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);
}
OpenPOWER on IntegriCloud