Merge from head@274682

4 files changed, 209 insertions, 10 deletions

diff --git a/games/primes/Makefile b/games/primes/Makefile
index 13c9048..bfc4147 100644
--- a/games/primes/Makefile
+++ b/games/primes/Makefile
@@ -2,7 +2,7 @@
 # $FreeBSD$
 
 PROG=	primes
-SRCS=	pattern.c pr_tbl.c primes.c
+SRCS=	pattern.c pr_tbl.c primes.c spsp.c
 MAN=
 DPADD=	${LIBM}
 LDADD=	-lm
diff --git a/games/primes/primes.c b/games/primes/primes.c
index 33838ee..a1c95c2 100644
--- a/games/primes/primes.c
+++ b/games/primes/primes.c
@@ -64,6 +64,7 @@ static const char rcsid[] =
 #include <ctype.h>
 #include <err.h>
 #include <errno.h>
+#include <inttypes.h>
 #include <limits.h>
 #include <math.h>
 #include <stdio.h>
@@ -111,10 +112,10 @@ main(int argc, char *argv[])
 	argv += optind;
 
 	start = 0;
-	stop = BIG;
+	stop = SPSPMAX;
 
 	/*
-	 * Convert low and high args.  Strtoul(3) sets errno to
+	 * Convert low and high args.  Strtoumax(3) sets errno to
 	 * ERANGE if the number is too large, but, if there's
 	 * a leading minus sign it returns the negation of the
 	 * result of the conversion, which we'd rather disallow.
@@ -126,18 +127,20 @@ main(int argc, char *argv[])
 			errx(1, "negative numbers aren't permitted.");
 
 		errno = 0;
-		start = strtoul(argv[0], &p, 0);
+		start = strtoumax(argv[0], &p, 0);
 		if (errno)
 			err(1, "%s", argv[0]);
 		if (*p != '\0')
 			errx(1, "%s: illegal numeric format.", argv[0]);
 
 		errno = 0;
-		stop = strtoul(argv[1], &p, 0);
+		stop = strtoumax(argv[1], &p, 0);
 		if (errno)
 			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. */
@@ -145,7 +148,7 @@ main(int argc, char *argv[])
 			errx(1, "negative numbers aren't permitted.");
 
 		errno = 0;
-		start = strtoul(argv[0], &p, 0);
+		start = strtoumax(argv[0], &p, 0);
 		if (errno)
 			err(1, "%s", argv[0]);
 		if (*p != '\0')
@@ -186,7 +189,7 @@ read_num_buf(void)
 		if (*p == '-')
 			errx(1, "negative numbers aren't permitted.");
 		errno = 0;
-		val = strtoul(buf, &p, 0);
+		val = strtoumax(buf, &p, 0);
 		if (errno)
 			err(1, "%s", buf);
 		if (*p != '\n')
@@ -241,7 +244,7 @@ primes(ubig start, ubig stop)
 		for (p = &prime[0], factor = prime[0];
 		    factor < stop && p <= pr_limit; factor = *(++p)) {
 			if (factor >= start) {
-				printf(hflag ? "0x%lx\n" : "%lu\n", factor);
+				printf(hflag ? "%" PRIx64 "\n" : "%" PRIu64 "\n", factor);
 			}
 		}
 		/* return early if we are done */
@@ -304,7 +307,11 @@ primes(ubig start, ubig stop)
 		 */
 		for (q = table; q < tab_lim; ++q, start+=2) {
 			if (*q) {
-				printf(hflag ? "0x%lx\n" : "%lu\n", start);
+				if (start > SIEVEMAX) {
+					if (!isprime(start))
+						continue;
+				}
+				printf(hflag ? "%" PRIx64 "\n" : "%" PRIu64 "\n", start);
 			}
 		}
 	}
diff --git a/games/primes/primes.h b/games/primes/primes.h
index 0951bbf..3a18fc7 100644
--- a/games/primes/primes.h
+++ b/games/primes/primes.h
@@ -41,8 +41,10 @@
  * chongo <for a good prime call: 391581 * 2^216193 - 1> /\oo/\
  */
 
+#include <stdint.h>
+
 /* ubig is the type that holds a large unsigned value */
-typedef unsigned long ubig;		/* must be >=32 bit unsigned value */
+typedef uint64_t ubig;			/* must be >=32 bit unsigned value */
 #define	BIG		ULONG_MAX	/* largest value will sieve */
 
 /* bytes in sieve table (must be > 3*5*7*11) */
@@ -57,6 +59,9 @@ typedef unsigned long ubig;		/* must be >=32 bit unsigned value */
 extern const ubig prime[];
 extern const ubig *const pr_limit;	/* largest prime in the prime array */
 
+/* Maximum size sieving alone can handle. */
+#define	SIEVEMAX 4295098368ULL
+
 /*
  * To avoid excessive sieves for small factors, we use the table below to
  * setup our sieve blocks.  Each element represents an odd number starting
@@ -64,3 +69,9 @@ extern const ubig *const pr_limit;	/* largest prime in the prime array */
  */
 extern const char pattern[];
 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/games/primes/spsp.c b/games/primes/spsp.c
new file mode 100644
index 0000000..f61acd6
--- /dev/null
+++ b/games/primes/spsp.c
@@ -0,0 +1,181 @@
+/*-
+ * Copyright (c) 2014 Colin Percival
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ *    notice, this list of conditions and the following disclaimer in the
+ *    documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+#include <sys/cdefs.h>
+__FBSDID("$FreeBSD$");
+
+#include <assert.h>
+#include <stddef.h>
+#include <stdint.h>
+
+#include "primes.h"
+
+/* Return a * b % n, where 0 <= a, b < 2^63, 0 < n < 2^63. */
+static uint64_t
+mulmod(uint64_t a, uint64_t b, uint64_t n)
+{
+	uint64_t x = 0;
+
+	while (b != 0) {
+		if (b & 1)
+			x = (x + a) % n;
+		a = (a + a) % n;
+		b >>= 1;
+	}
+
+	return (x);
+}
+
+/* Return a^r % n, where 0 <= a < 2^63, 0 < n < 2^63. */
+static uint64_t
+powmod(uint64_t a, uint64_t r, uint64_t n)
+{
+	uint64_t x = 1;
+
+	while (r != 0) {
+		if (r & 1)
+			x = mulmod(a, x, n);
+		a = mulmod(a, a, n);
+		r >>= 1;
+	}
+
+	return (x);
+}
+
+/* Return non-zero if n is a strong pseudoprime to base p. */
+static int
+spsp(uint64_t n, uint64_t p)
+{
+	uint64_t x;
+	uint64_t r = n - 1;
+	int k = 0;
+
+	/* Compute n - 1 = 2^k * r. */
+	while ((r & 1) == 0) {
+		k++;
+		r >>= 1;
+	}
+
+	/* Compute x = p^r mod n.  If x = 1, n is a p-spsp. */
+	x = powmod(p, r, n);
+	if (x == 1)
+		return (1);
+
+	/* Compute x^(2^i) for 0 <= i < n.  If any are -1, n is a p-spsp. */
+	while (k > 0) {
+		if (x == n - 1)
+			return (1);
+		x = powmod(x, 2, n);
+		k--;
+	}
+
+	/* Not a p-spsp. */
+	return (0);
+}
+
+/* Test for primality using strong pseudoprime tests. */
+int
+isprime(ubig _n)
+{
+	uint64_t n = _n;
+
+	/*
+	 * Values from:
+	 * C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.,
+	 * The pseudoprimes to 25 * 10^9, Math. Comp. 35(151):1003-1026, 1980.
+	 */
+
+	/* No SPSPs to base 2 less than 2047. */
+	if (!spsp(n, 2))
+		return (0);
+	if (n < 2047ULL)
+		return (1);
+
+	/* No SPSPs to bases 2,3 less than 1373653. */
+	if (!spsp(n, 3))
+		return (0);
+	if (n < 1373653ULL)
+		return (1);
+
+	/* No SPSPs to bases 2,3,5 less than 25326001. */
+	if (!spsp(n, 5))
+		return (0);
+	if (n < 25326001ULL)
+		return (1);
+
+	/* No SPSPs to bases 2,3,5,7 less than 3215031751. */
+	if (!spsp(n, 7))
+		return (0);
+	if (n < 3215031751ULL)
+		return (1);
+
+	/*
+	 * Values from:
+	 * G. Jaeschke, On strong pseudoprimes to several bases,
+	 * Math. Comp. 61(204):915-926, 1993.
+	 */
+
+	/* No SPSPs to bases 2,3,5,7,11 less than 2152302898747. */
+	if (!spsp(n, 11))
+		return (0);
+	if (n < 2152302898747ULL)
+		return (1);
+
+	/* No SPSPs to bases 2,3,5,7,11,13 less than 3474749660383. */
+	if (!spsp(n, 13))
+		return (0);
+	if (n < 3474749660383ULL)
+		return (1);
+
+	/* No SPSPs to bases 2,3,5,7,11,13,17 less than 341550071728321. */
+	if (!spsp(n, 17))
+		return (0);
+	if (n < 341550071728321ULL)
+		return (1);
+
+	/* No SPSPs to bases 2,3,5,7,11,13,17,19 less than 341550071728321. */
+	if (!spsp(n, 19))
+		return (0);
+	if (n < 341550071728321ULL)
+		return (1);
+
+	/*
+	 * Value from:
+	 * Y. Jiang and Y. Deng, Strong pseudoprimes to the first eight prime
+	 * bases, Math. Comp. 83(290):2915-2924, 2014.
+	 */
+
+	/* No SPSPs to bases 2..23 less than 3825123056546413051. */
+	if (!spsp(n, 23))
+		return (0);
+	if (n < 3825123056546413051)
+		return (1);
+
+	/* We can't handle values larger than this. */
+	assert(n <= SPSPMAX);
+
+	/* UNREACHABLE */
+	return (0);
+}