svn commit: r320218 - in stable/11/usr.bin: factor primes

[Colin Percival]

Author: cperciva
Date: Thu Jun 22 05:26:08 2017
New Revision: 320218
URL: https://svnweb.freebsd.org/changeset/base/320218

Log:
  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

Modified:
  stable/11/usr.bin/factor/factor.6
  stable/11/usr.bin/primes/primes.c
  stable/11/usr.bin/primes/primes.h
  stable/11/usr.bin/primes/spsp.c
Directory Properties:
  stable/11/   (props changed)

Modified: stable/11/usr.bin/factor/factor.6
==============================================================================
--- stable/11/usr.bin/factor/factor.6	Thu Jun 22 05:22:21 2017	(r320217)
+++ stable/11/usr.bin/factor/factor.6	Thu Jun 22 05:26:08 2017	(r320218)
@@ -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.

Modified: stable/11/usr.bin/primes/primes.c
==============================================================================
--- stable/11/usr.bin/primes/primes.c	Thu Jun 22 05:22:21 2017	(r320217)
+++ stable/11/usr.bin/primes/primes.c	Thu Jun 22 05:26:08 2017	(r320218)
@@ -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. */

Modified: stable/11/usr.bin/primes/primes.h
==============================================================================
--- stable/11/usr.bin/primes/primes.h	Thu Jun 22 05:22:21 2017	(r320217)
+++ stable/11/usr.bin/primes/primes.h	Thu Jun 22 05:26:08 2017	(r320218)
@@ -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: stable/11/usr.bin/primes/spsp.c
==============================================================================
--- stable/11/usr.bin/primes/spsp.c	Thu Jun 22 05:22:21 2017	(r320217)
+++ stable/11/usr.bin/primes/spsp.c	Thu Jun 22 05:26:08 2017	(r320218)
@@ -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);
 }