Imported and adapted from FreeBSD svn r272166 and r272207; this fixes

false positives for products of primes larger than 2^16. For example,
before this commit:

  $ /usr/games/primes 4295360521 4295360522
  4295360521
but
  $ /usr/games/factor 4295360521
  4295360521: 65539 65539

or
  $ /usr/games/primes 3825123056546413049 3825123056546413050
  3825123056546413049
yet
  $ /usr/games/factor 3825123056546413049
  3825123056546413049: 165479 23115459100831

or
  $ /usr/games/primes 18446744073709551577
  18446744073709551577
although
  $ /usr/games/factor 18446744073709551577
  18446744073709551577: 139646831 132095686967

Incidentally, the above examples show the smallest and largest cases that
were erroneously stated as prime in the range 2^32 .. 3825123056546413049
.. 2^64; the primes(6) program now stops at 3825123056546413050 as
primality tests on larger integers would be by brute force factorization.

In addition, special to the NetBSD version:
. for -d option, skip first difference when start is >65537 as it is incorrect
. corrected usage to mention both the existing -d as well as the new -h option

For original FreeBSD commit message by Colin Percival, see:
http://svnweb.freebsd.org/base?view=revision&revision=272166

1 files changed, 195 insertions, 0 deletions

diff --git a/primes/spsp.c b/primes/spsp.c
new file mode 100644
index 00000000..7d1b0841
--- /dev/null
+++ b/primes/spsp.c
@@ -0,0 +1,195 @@
+/*	$NetBSD: spsp.c,v 1.1 2014/10/02 21:36:37 ast Exp $	*/
+
+/*-
+ * 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>
+#ifndef lint
+__COPYRIGHT("@(#) Copyright (c) 1989, 1993\
+ The Regents of the University of California.  All rights reserved.");
+#endif /* not lint */
+
+#ifndef lint
+#if 0
+static char sccsid[] = "@(#)primes.c    8.5 (Berkeley) 5/10/95";
+#else
+__RCSID("$NetBSD: spsp.c,v 1.1 2014/10/02 21:36:37 ast Exp $");
+#endif
+#endif /* not lint */
+
+#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(uint64_t _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);
+}