[src/trunk]: src/games/primes PR/52976: Eitan Adler: handle larger primes

[christos <christos%NetBSD.org@localhost>]

details:   https://anonhg.NetBSD.org/src/rev/d0cdf33b2f86
branches:  trunk
changeset: 829478:d0cdf33b2f86
user:      christos <christos%NetBSD.org@localhost>
date:      Sat Feb 03 15:40:29 2018 +0000

description:
PR/52976: Eitan Adler: handle larger primes
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

diffstat:

 games/primes/primes.6 |  13 +++----------
 games/primes/primes.c |   9 +++------
 games/primes/primes.h |   5 +----
 games/primes/spsp.c   |  44 ++++++++++++++++++++++++++++++++------------
 4 files changed, 39 insertions(+), 32 deletions(-)

diffs (170 lines):

diff -r 9816aecea86a -r d0cdf33b2f86 games/primes/primes.6
--- a/games/primes/primes.6     Sat Feb 03 11:30:01 2018 +0000
+++ b/games/primes/primes.6     Sat Feb 03 15:40:29 2018 +0000
@@ -1,4 +1,4 @@
-.\"    $NetBSD: primes.6,v 1.5 2014/10/04 13:15:50 wiz Exp $
+.\"    $NetBSD: primes.6,v 1.6 2018/02/03 15:40:29 christos Exp $
 .\"
 .\" Copyright (c) 1989, 1993
 .\"    The Regents of the University of California.  All rights reserved.
@@ -35,7 +35,7 @@
 .\"
 .\" By Landon Curt Noll, http://www.isthe.com/chongo/index.html /\oo/\
 .\"
-.Dd February 3, 2008
+.Dd February 2, 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.
diff -r 9816aecea86a -r d0cdf33b2f86 games/primes/primes.c
--- a/games/primes/primes.c     Sat Feb 03 11:30:01 2018 +0000
+++ b/games/primes/primes.c     Sat Feb 03 15:40:29 2018 +0000
@@ -1,4 +1,4 @@
-/*     $NetBSD: primes.c,v 1.21 2014/10/04 13:15:50 wiz Exp $  */
+/*     $NetBSD: primes.c,v 1.22 2018/02/03 15:40:29 christos Exp $     */
 
 /*
  * Copyright (c) 1989, 1993
@@ -42,7 +42,7 @@
 #if 0
 static char sccsid[] = "@(#)primes.c   8.5 (Berkeley) 5/10/95";
 #else
-__RCSID("$NetBSD: primes.c,v 1.21 2014/10/04 13:15:50 wiz Exp $");
+__RCSID("$NetBSD: primes.c,v 1.22 2018/02/03 15:40:29 christos Exp $");
 #endif
 #endif /* not lint */
 
@@ -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. */
diff -r 9816aecea86a -r d0cdf33b2f86 games/primes/primes.h
--- a/games/primes/primes.h     Sat Feb 03 11:30:01 2018 +0000
+++ b/games/primes/primes.h     Sat Feb 03 15:40:29 2018 +0000
@@ -1,4 +1,4 @@
-/*     $NetBSD: primes.h,v 1.6 2014/10/02 21:36:37 ast Exp $   */
+/*     $NetBSD: primes.h,v 1.7 2018/02/03 15:40:29 christos Exp $      */
 
 /*
  * Copyright (c) 1989, 1993
@@ -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
diff -r 9816aecea86a -r d0cdf33b2f86 games/primes/spsp.c
--- a/games/primes/spsp.c       Sat Feb 03 11:30:01 2018 +0000
+++ b/games/primes/spsp.c       Sat Feb 03 15:40:29 2018 +0000
@@ -1,4 +1,4 @@
-/*     $NetBSD: spsp.c,v 1.1 2014/10/02 21:36:37 ast Exp $     */
+/*     $NetBSD: spsp.c,v 1.2 2018/02/03 15:40:29 christos Exp $        */
 
 /*-
  * Copyright (c) 2014 Colin Percival
@@ -36,7 +36,7 @@
 #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 $");
+__RCSID("$NetBSD: spsp.c,v 1.2 2018/02/03 15:40:29 christos Exp $");
 #endif
 #endif /* not lint */
 
@@ -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);
-
-       /* 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);
 }