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
This commit is contained in:
christos
parent
317df2f20a
commit
2429b427fa
@ -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 @@ Originally by
|
||||
.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.
|
||||
|
||||
@ -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 @@ __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
|
||||
#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 @@ 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
|
||||
@ -145,9 +145,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 (>%" PRIu64 ").",
|
||||
argv[1], (uint64_t) SPSPMAX);
|
||||
break;
|
||||
case 1:
|
||||
/* Start on the command line. */
|
||||
|
||||
@ -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 @@ extern const size_t pattern_size; /* length of pattern array */
|
||||
|
||||
/* Test for primality using strong pseudoprime tests. */
|
||||
int isprime(uint64_t);
|
||||
|
||||
/* Maximum value which the SPSP code can handle. */
|
||||
#define SPSPMAX 3825123056546413050ULL
|
||||
|
||||
@ -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 @@ __COPYRIGHT("@(#) Copyright (c) 1989, 1993\
|
||||
#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 @@ __RCSID("$NetBSD: spsp.c,v 1.1 2014/10/02 21:36:37 ast Exp $");
|
||||
|
||||
#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 @@ isprime(uint64_t _n)
|
||||
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);
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user