Correctly enumerate primes between 4295098369 and 3825123056546413050.
Prior to this commit, primes(6) relied solely on sieving with primes up to 65537, with the effect that composite numbers which are the product of two non-16-bit primes would be incorrectly identified as prime. For example, # primes 1099511627800 1099511627820 would output 1099511627803 1099511627807 1099511627813 when in fact only the first of those values is prime. This commit adds strong pseudoprime tests to validate the candidates which pass the initial sieving stage, using bases of 2, 3, 5, 7, 11, 13, 17, 19, and 23. Thanks to papers from C. Pomerance, J.L. Selfridge, and S.S. Wagstaff, Jr.; G. Jaeschke; and Y. Jiang and Y. Deng, we know that the smallest value which passes these tests is 3825123056546413051. At present we do not know how many strong pseudoprime tests are required to prove primality for values larger than 3825123056546413050, so we force primes(6) to stop at that point. Reviewed by: jmg Relnotes: primes(6) now correctly enumerates primes up to 3825123056546413050 MFC after: 7 days Sponsored by: EuroBSDCon devsummit
This commit is contained in:
parent
5067548af0
commit
535ab8fde8
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=272166
@ -90,7 +90,7 @@ value must not be greater than the maximum.
|
||||
The default and maximum value of
|
||||
.Ar stop
|
||||
is 4294967295 on 32-bit architectures
|
||||
and 18446744073709551615 on 64-bit ones.
|
||||
and 3825123056546413050 on 64-bit ones.
|
||||
.Pp
|
||||
When the
|
||||
.Nm primes
|
||||
@ -120,3 +120,9 @@ cannot handle the
|
||||
factor list,
|
||||
.Nm primes
|
||||
will not get you a world record.
|
||||
.Pp
|
||||
.Nm primes
|
||||
is unable to list primes between 3825123056546413050 and 18446744073709551615
|
||||
since it relies on strong pseudoprime tests after sieving, and nobody has
|
||||
proven how many strong pseudoprime tests are required to prove primality for
|
||||
integers larger than 3825123056546413050.
|
||||
|
@ -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
|
||||
|
@ -111,7 +111,7 @@ main(int argc, char *argv[])
|
||||
argv += optind;
|
||||
|
||||
start = 0;
|
||||
stop = BIG;
|
||||
stop = (sizeof(ubig) > 4) ? SPSPMAX : BIG;
|
||||
|
||||
/*
|
||||
* Convert low and high args. Strtoul(3) sets errno to
|
||||
@ -138,6 +138,8 @@ main(int argc, char *argv[])
|
||||
err(1, "%s", argv[1]);
|
||||
if (*p != '\0')
|
||||
errx(1, "%s: illegal numeric format.", argv[1]);
|
||||
if ((uint64_t)stop > SPSPMAX)
|
||||
errx(1, "%s: stop value too large.", argv[1]);
|
||||
break;
|
||||
case 1:
|
||||
/* Start on the command line. */
|
||||
@ -304,6 +306,10 @@ primes(ubig start, ubig stop)
|
||||
*/
|
||||
for (q = table; q < tab_lim; ++q, start+=2) {
|
||||
if (*q) {
|
||||
if ((uint64_t)start > SIEVEMAX) {
|
||||
if (!isprime(start))
|
||||
continue;
|
||||
}
|
||||
printf(hflag ? "0x%lx\n" : "%lu\n", start);
|
||||
}
|
||||
}
|
||||
|
@ -57,6 +57,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 +67,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
|
||||
|
181
games/primes/spsp.c
Normal file
181
games/primes/spsp.c
Normal file
@ -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);
|
||||
}
|
Loading…
Reference in New Issue
Block a user