bin/52976: [primes] handle larger primes

[lists%eitanadler.com@localhost]

>Number:         52976
>Category:       bin
>Synopsis:       [primes] handle larger primes
>Confidential:   no
>Severity:       non-critical
>Priority:       low
>Responsible:    bin-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Sat Feb 03 03:40:00 +0000 2018
>Originator:     Eitan Adler
>Release:        HEAD
>Organization:
>Environment:
>Description:
Prime numbers are cool. We'd like more of them.

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.

Adapted from: FreeBSD
>How-To-Repeat:

>Fix:
? a.out
? obj
Index: primes.6
===================================================================
RCS file: /cvsroot/src/games/primes/primes.6,v
retrieving revision 1.5
diff -u -r1.5 primes.6
--- primes.6	4 Oct 2014 13:15:50 -0000	1.5
+++ primes.6	3 Feb 2018 03:30:29 -0000
@@ -35,7 +35,7 @@
 .\"
 .\" By Landon Curt Noll, http://www.isthe.com/chongo/index.html /\oo/\
 .\"
-.Dd February 3, 2008
+.Dd February 02, 2018
 .Dt PRIMES 6
 .Os
 .Sh NAME
@@ -100,14 +100,7 @@
 .An Landon Curt Noll ,
 extended to some 64-bit primes by
 .An Colin Percival .
-.Sh CAVEATS
+.Sh BUGS
 This
 .Nm
 program won't get you a world record.
-.Pp
-The program is not able to list primes between
-3825123056546413050 and 18446744073709551615 (2^64
-- 1) as it relies on strong pseudoprime tests after
-sieving, and it is yet unknown how many of those
-tests are needed to prove primality for integers
-larger than 3825123056546413050.
Index: primes.c
===================================================================
RCS file: /cvsroot/src/games/primes/primes.c,v
retrieving revision 1.21
diff -u -r1.21 primes.c
--- primes.c	4 Oct 2014 13:15:50 -0000	1.21
+++ primes.c	3 Feb 2018 03:30:29 -0000
@@ -118,7 +118,7 @@
 	argv += optind;
 
 	start = 0;
-	stop = SPSPMAX;
+	stop = (uint64_t)(-1);
 
 	/*
 	 * Convert low and high args.  Strtoumax(3) sets errno to
@@ -145,9 +145,6 @@
 			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 (>%" PRIu64 ").",
-				argv[1], (uint64_t) SPSPMAX);
 		break;
 	case 1:
 		/* Start on the command line. */
Index: primes.h
===================================================================
RCS file: /cvsroot/src/games/primes/primes.h,v
retrieving revision 1.6
diff -u -r1.6 primes.h
--- primes.h	2 Oct 2014 21:36:37 -0000	1.6
+++ primes.h	3 Feb 2018 03:30:29 -0000
@@ -69,6 +69,3 @@
 
 /* Test for primality using strong pseudoprime tests. */
 int isprime(uint64_t);
-
-/* Maximum value which the SPSP code can handle. */
-#define	SPSPMAX 3825123056546413050ULL
Index: spsp.c
===================================================================
RCS file: /cvsroot/src/games/primes/spsp.c,v
retrieving revision 1.1
diff -u -r1.1 spsp.c
--- spsp.c	2 Oct 2014 21:36:37 -0000	1.1
+++ spsp.c	3 Feb 2018 03:30:29 -0000
@@ -46,23 +46,33 @@
 
 #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)
 {
@@ -186,10 +196,20 @@
 		return (0);
 	if (n < 3825123056546413051)
 		return (1);
+	/*
+	 * Value from:
+	 * J. Sorenson and J. Webster, Strong pseudoprimes to twelve prime
+	 * bases, Math. Comp. 86(304):985-1003, 2017.
+	 */
 
-	/* We can't handle values larger than this. */
-	assert(n <= SPSPMAX);
+       /* 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);
 
-	/* UNREACHABLE */
-	return (0);
+       /* All 64-bit values are less than 318665857834031151167461. */
+       return (1);
 }