Let be an integer, and let …, , be admissible linear polynomials over the integers, or the pattern. We present two algorithms that find all integers where and all the are prime.
Our first algorithm takes arithmetic operations using space.
Our second algorithm takes slightly more time, arithmetic operations, but uses only space for a fixed constant. The correctness of this algorithm is unconditional, but our analysis of its running time depends on two reasonable but unproven conjectures.
We are unaware of any previous complexity results for this problem beyond the use of a prime sieve.
We also implemented several parallel versions of our second algorithm to show it is viable in practice. In particular, we found some new Cunningham chains of length 15, and we found all quadruplet primes up to .
Mathematicians have long been interested in prime numbers and how they appear in patterns. (See, for example, Reference 12, Ch. A.) In this paper, we are interested in the complexity of the following algorithmic problem:
Given a pattern and a bound find all primes , that fit the pattern.
To address this, first we will discuss and define a pattern of primes, then we will look at what is known about the distribution of primes in patterns to see what we can reasonably expect for the complexity of this problem, and finally we will discuss previous work and state our new results.
1.1. Prime patterns
Perhaps the simplest of patterns of primes are the twin primes, which satisfy the pattern where both and are prime. Examples include 59, 61 and 101, 103.
We can, of course, generalize this to larger patterns. For example, prime quadruplets have the form , , , and examples include 11, 13, 17, 19 and 1481, 1483, 1487, 1489. ,
Larger patterns of primes of this type are called prime . If the -tuples has the smallest possible difference between its first and last primes (its diameter), it is also called a prime constellation. The pattern -tuple , , , is a constellation, as is the smallest possible diameter for a pattern of length 4. There are two constellations of length 3: , , and , , See, for example, .Reference 7, §1.2.2 or Reference 25, Ch. 3.
Sophie Germain studied the pattern which was later generalized to Cunningham chains of two kinds. Chains of the first kind have the pattern , , , , , and chains of the second kind have the pattern , , , , ,.
Chernick Reference 6 showed that any prime pattern of the form , , gives a Carmichael number composed of the product of these three primes.
Let be an integer. A prime pattern of size is a list of linear polynomials over the integers with positive leading coefficients, A pattern of size . is admissible if for every prime there is an integer , such that does not divide any of the For an algorithm to test for pattern admissibility, see .Reference 25, pp. 62–63.
We restrict our notion of pattern to linear polynomials in this paper.
1.2. The distribution of primes in patterns
The unproven twin prime conjecture states there are infinitely many twin primes. Yitang Zhang Reference 33 recently showed that there is a positive integer such that the pattern is satisfied by infinitely many primes.
The Hardy-Littlewood conjecture -tupleReference 14 implies that each pattern, with leading coefficients of 1, that is admissible will be satisfied by primes infinitely often. Further, the conjecture implies that the number of primes in such a pattern of length is roughly proportional to .
For twin primes, then, the Hardy-Littlewood conjecture gives an estimate of
for the number of twin primes where , is the twin primes constant. Brun’s theorem gives an upper bound for the number of twin primes For every . with -tuple there is a corresponding conjectured constant of proportionality that depends on the pattern and on , generalizing the twin primes constant. ,
Dickson’s conjecture Reference 8 states that there are infinitely many primes satisfying any fixed admissible pattern, even with leading coefficients Thus, it applies to fixed-length Cunningham chains as well. .
We have the following upper bound, due to Halberstam and Richert Reference 13, Thm. 2.4.
Given an admissible pattern of length the number of integers , such that the are all simultaneously prime and is .
Here the implied constant of the big- can depend on but does not depend on the coefficients of the .
1.3. Previous work
We can obtain a rather straightforward analysis by simply finding all primes and then scanning to see how many tuples match the pattern. Note that it is not necessary to write all the primes down; we can scan them in small batches as we go. Since the scan phase is fast, the bottleneck would be finding the primes using a sieve. The Atkin-Bernstein sieve Reference 2 does this using at most arithmetic operations and space. Galway Reference 10 showed how to reduce space further to roughly but this version could not use the wheel data structure and requires , arithmetic operations. See also Reference 16.
We follow the convention from the analysis of prime sieves of not charging for the space used by the output of the algorithm. Note that if we charged for the space of the output, we could not expect to do better than bits in general.
A more specialized algorithm that searches for only the primes in the pattern and not all the primes should do much better. Günter Löh Reference 20 and Tony Forbes Reference 9 described algorithms to do exactly this but gave no runtime analysis. It seems likely their algorithms are faster than but without an analysis, we don’t know if this is true or not. Note that Forbes outlined an odometer-like data structure that seems to be similar to the wheel data structure we employ. ,
Of course, by the prime number theorem, there are only about primes so the current best prime sieves use , arithmetic operations per prime on average. We do not know if anything smaller is possible. Applying this average cost to the results of Lemma 1, we can hope for an algorithm that takes arithmetic operations to find all primes in a fixed pattern of length .
In essence, this is what we prove as our main result.
1.4. New results
Our contribution is the following.
There is an algorithm that when given a list of distinct linear polynomials over the integers, with positive leading coefficients, (the pattern), and a search bound finds all integers such that and all the are prime. This algorithm uses at most arithmetic operations and bits of space.
The space needed by this algorithm limits its practicality. By replacing the Atkin-Bernstein sieve with the sieve of Eratosthenes combined with prime tests, we can greatly reduce the need for space.
Let There is an algorithm that when given a list of . distinct linear polynomials over the integers, with positive leading coefficients, (the pattern), and a search bound finds all integers such that and all the are prime. This algorithm uses at most arithmetic operations and bits of space, under the assumption of Conjectures 1 and 2 (see below). Correctness of the output is unconditional.
If then Conjecture ,1 is not needed. We use the sieve of Eratosthenes with primes up to a bound after which we apply a base-2 pseudoprime test and then a version of the AKS prime test ,Reference 1 that was improved to take time Reference 19. Our goal here is to balance the cost of sieving with the cost of prime testing. For the range to keep the cost of prime testing low enough, we replace the AKS prime test with the pseudosquares prime test of Lukes, Patterson, and Williams ,Reference 21. This prime test takes only time under the condition of an unproven but reasonable conjecture on the distribution of pseudosquares due to Bach and Huelsbergen Reference 3. Note that the correctness of our algorithm does not rely on any unproven conjectures.
The pseudosquare is the smallest positive integer that is not a square, is and is a quadratic residue modulo every odd prime ,.
Conjecture 1 (Bach and Huelsbergen Reference 3).
Let be the largest pseudosquare Then ..
Our second conjecture is needed to bound the cost of prime testing after sieving out by primes in an arithmetic progression. Let denote the smallest prime divisor of the integer .
Let be positive integers with Then .
We performed a few computations with this second version of our algorithm to show its practicality. A couple of these computational results are new.
The rest of this paper is organized as follows. In §2 we present our proof of Theorem 1, including our model of computation in §2.1, a description of our first algorithm in §2.2, and its running time analysis in §2.3. In §3 we discuss our second algorithm in §3.1 and its analysis in §3.2, thereby proving Theorem 2. We present our computational results in §4, including our work on twin primes in §4.1, our work on prime quadruplets in §4.2, and our results on Cunningham chains in §4.3. We wrap up with a discussion of possible future work in §5.
2.1. Model of computation
Our model of computation is a standard random access machine with infinite, direct-access memory. Memory can be addressed at the bit level or at the word level, and the word size is bits if is the input. Arithmetic operations on integers of bits take constant time, as do memory/array accesses, comparisons, and other basic operations.
We count space used in bits, and we do not include the size of the output.
2.2. Our first algorithm
In this section we present the version of our algorithm with the smallest running time; we perform the analysis in the next section.
The input to the algorithm is the search bound and the pattern, which consists of the value of and the list of linear polynomials We write . for the multiplier and for the offset for each form For simplicity, we often assume that . and but this convenience is not required to obtain our complexity bound. So, for example, for Cunningham chains of the first kind we would have , and .
We begin by finding the list of primes up to and choosing a subset to be the wheel primes Define . Generally, we put all primes up to a bound . into with , maximal such that Then by the prime number theorem .Reference 15 we have
This implies that .
In practice, if there is a prime that provides poor filtering for the pattern, we consider dropping it from and increasing to include another prime.
We must also check all primes to see if they participate in the prime pattern.
Next, we construct the wheel data structure so that it will generate all acceptable residues modulo .
The data structure is initialized with a list of pairwise coprime moduli and for each such modulus a list of acceptable residues mod , encoded as a bit vector of length , The wheel modulus . is then the product of these moduli. Once initialized, the wheel has the ability to enumerate suitable residues modulo in amortized constant time per residue. The residues do not appear in any particular order.
Therefore, for each prime we compute a bit vector (ones) that encodes the list of acceptable residues. For any integral we want , to not be divisible by So if . divides or, equivalently, if does not divide and so then bit position , must be a zero.
Continuing the example above, for Cunningham chains of the first kind, gives the vector 001, and gives the vector 0010111.
We then construct the wheel data structure as described in Reference 28, §4.
For each residue generated by the wheel, we sieve arithmetic progressions for primes up to , or , for …, , We do this using the Atkin-Bernstein sieve. (See .Reference 2, §5 for how to use the sieve to find primes in an arithmetic progression.)
For each this yields , bit vectors of length which are combined using a bitwise AND operation to obtain the bit positions for where the pattern is satisfied by primes.
Let us find the prime quadruplets ( , up to ) using a wheel of size .
We generate just the one quadruplet that contains wheel primes.
The only acceptable residue mod is for , it is (if then , is divisible by giving ), and for , it is giving , For . there are three acceptable residues, giving us , by the Chinese remainder theorem.
For we sieve the progression , for primes, getting the bit vector 10110; 11, 431, and 641 are prime, 221 and 851 are not. We then sieve for primes, getting 11111; they are all prime. Sieving gives 11011, and gives 11101.
Bitwise AND-ing these four vectors together gives 10000. The only quadruplet we find is , , ,.
Doing the same for we get the following: ,
which gives us the quadruplet , , ,.
Finally for we get
and we get two quadruplets, , , , and , , ,.
2.3. Complexity analysis
We’ll look at the cost for each step of the algorithm above.
We can use the Atkin-Bernstein sieve to find the primes up to in arithmetic operations using space.
Recall that the largest wheel prime is roughly Constructing the bit vector .ones for one prime takes time to initially write all ones and then time to mark the zeros. Summing over all wheel primes gives operations.
From Reference 28, Thm. 4.1 the total cost to build the wheel is operations, and it occupies space.
The Atkin-Bernstein sieve finds all primes in an arithmetic progression in an interval of size or larger in time linear in the length of the interval using space proportional to Reference 2, §5. Therefore, sieving for primes takes operations for each of the residue classes for a total of , The cost to generate each value of . using the wheel is negligible in comparison. The space used is bits.
Next we show that the total number of residues is roughly asymptotic to For a pattern . of size all but finitely many primes , will have possible residues. Let be a constant, depending on the pattern such that all primes larger than , have residues. Recall that is a bound for the largest prime in the wheel. Then the total number of residues will be asymptotically bounded by
By Bernoulli’s inequality we have so that ,
The first product is bounded by a constant that depends on which we refer to as , going forward. If then by ,Reference 26, (3.26) we know that
As shown above, if is the largest prime in we have , Asymptotically, . is a very safe underestimate. Pulling this all together, we obtain the bound
for the total number of residues .
Multiplying the sieving cost by the number of residues gives
We have proved Theorem 1.
It is possible to pin down the constant that depends on , (and We can use ).Reference 26, (3.30) and, assuming (note that we have ),
Bringing in the factor of introduced in the sieving, we have
It remains to get a bound for Set . an upper bound on all the coefficients in the pattern. Set , For a prime . if , has repeated roots modulo then , But then . for some with which means , divides and hence , Thus . and we have the bound ,
For specific patterns, the constant can be computed explicitly, typically giving much better results than this. In particular, we assumed every prime contributed all its residues (which happens if a prime divides all the which is unusual. ’s),
As an example, for our pattern we have , and our constant would be ,
Here the and multipliers are included because and have only one residue each instead of and respectively. ,
Note that one of the anonymous referees deserves the primary credit for this bound on .
The primary difficulty in reaching large values of with our first algorithm is the amount of space it requires. One way to address this is to create a larger wheel, sieve more but shorter arithmetic progressions for primes, and rely less on sieving and more on primality tests (in the style of Reference 28) when searching for -tuples.
We use the sieve of Eratosthenes instead of the Atkin-Bernstein sieve for the arithmetic progressions, and this is the source of the factor slowdown. The gains here are less space needed by a factor of and the effective trial division performs a quantifiable filtering out of non-primes. ,
Instead of sieving by all primes up to we sieve only by primes up to a bound , for a constant In practice, we choose . so everything, including space for the wheel data structure, the bit vector for sieving, and the list of primes fits in CPU cache. We then use a base-2 pseudoprime test, followed by a prime test as needed. For smaller , we use the pseudosquares prime test of Lukes, Patterson, and Williams ,Reference 21, which is fast and deterministic, assuming a sufficient table of pseudosquares is available. Importantly, it takes advantage of the trial division effect of the sieve of Eratosthenes. For larger we can simply use the AKS prime test ,Reference 1.
This change means we can get by with only space. Choosing larger or smaller gives a tradeoff between the cost of sieving and the cost of performing base-2 pseudoprime tests.
3.1. Our second algorithm
Choose a constant and then set , a power of , near We begin by finding the list of primes up to . and dividing them into the two sets and Small primes go into . and the remainder go in We want . to be as large as possible with the constraint that This implies that the largest prime in . will be roughly .
If we will need to perform the pseudosquares prime test, so in preparation, find all pseudosquares , (see Reference 21).
Next, as before, we construct the wheel data structure so that it will generate all possible correct residues modulo .
For each residue generated by the wheel, we construct a bit vector v of length Each vector position .v, for …, , represents , for the -tuple , …, , We initialize .v but clear it to if we find a prime where for some .
Once this sieving is complete, the only integers with v that remain satisfy the property that all the have no prime divisors less than .
For each such remaining (that is, v we first do a base-2 strong pseudoprime test on ), If it fails, we cross it off (set .v If it passes, we try ). and so forth, keeping v only if all values pass the pseudoprime test. We then perform a full prime test on the for all If . we use the Lukes, Patterson, and Williams pseudosquares prime test ,Reference 21 as done in Reference 28. For larger we use the AKS prime test ,Reference 1. (This is for the purposes of the theorem; in practice, the pseudosquares prime test is faster, so we use that instead.) If all the pass the prime tests, the corresponding is written for output. -tuple
This version of the algorithm works best for When . the prime tests become the runtime bottleneck, and so we recommend using , so that the base-2 pseudoprime tests and the pseudosquares prime test are not needed, as the sieving will leave only primes.
We use the same example as above, finding prime quadruplets up to which uses the pattern , , , , We go a bit higher this time to illustrate the sieving. Recall that our wheel uses the primes . and generates the three residues , , modulo .
We will use as our sieve bound, so , , , As with the wheel primes, quadruplets that include sieve primes must be generated separately, so we output the quadruplet . , , , at this point.
For we have our progression ,
For we cross off integers that are , , , , modulo This takes four passes: on the first pass we remove . , , which are , on the second pass we remove ; , which are , (for example the third removes ); , which are , and the fourth removes ; , which are , Observe that for a given residue modulo . the numbers to cross off are exactly , spaces apart.
For we cross off integers that are , , , , modulo We remove . , for , , for , for and , , for .
For we cross off integers that are , , , , modulo We remove . , for , for , for and , , for .
For we cross off integers that are , , , , modulo We remove . , for , , for , for and , for .
At this point we perform base-2 strong pseudoprime tests, followed by prime tests as needed. Here 851 and 3161 are not prime, and 1481 leads to the quadruplet , , ,.
This is then repeated with ,.
3.2. Complexity analysis
Finding the primes up to takes time using space, well within our bounds. See Reference 28 for a sublinear time algorithm to find all needed pseudosquares. In practice, all pseudosquares up to are known Reference 29. The cost in time and space to build the wheel is, up to a constant factor, the same. So we now focus on steps (4) and (5).
As shown above, the number of residues to check is The time to sieve each interval of length . using primes up to is at most proportional to
Here the multiplier is required because we cross off residues modulo most of the primes That said, this multiple of . can be absorbed into the implied constant that depends on the pattern , from earlier. ,
At this point we make use of Conjecture 2 to bound the number of integers free from prime divisors in an arithmetic progression.
With this assumption in hand, by Mertens’s theorem, an average of at most
vector locations remain to be prime tested. (Note that we cannot make any assumptions about the relative independence of the primality of the values for different and so we cannot use a , factor here.)
A single base-2 strong pseudoprime test takes operations to perform, for a total cost proportional to
arithmetic operations to do the base-2 strong pseudoprime tests for each value of This matches the sieving cost of . from above. (Note that if we deliberately choose a larger value for the increased sieving will decrease the number of pseudoprime tests needed. This tradeoff can be used to fine-tune the running time of the algorithm.) ,
Thus, the total cost for sieving and base-2 pseudoprime tests is
which we obtain by multiplying by the number of residues .
Next we need to count integers that pass the base-2 strong pseudoprime test. Such integers are either prime or composite base-2 pseudoprimes. We switch to counting across all residues to obtain an overall bound.
Lemma 1 tells us that at most integers are prime that fit the pattern, so this is an upper bound on primes that pass the base-2 pseudoprime test.
Pomerance Reference 24 showed that the number of composite base-2 pseudoprimes is bounded by
which is negligible. This plus the bound for primes above gives us the bound we desire for all integers that pass the base-2 pseudoprime test.
Next, to bound the cost of prime tests, we have two cases: or .
For we use the AKS prime test ,Reference 1 (improved in Reference 19) which takes time The cost of applying the AKS prime test to all the integers . after they all pass a base-2 pseudoprime test is at most proportional to
which is for .
Note that when is large, in practice we might only do the base-2 pseudoprime tests and then run full prime tests on the output afterwards, since the amount of output will be rather small.
For Conjecture ,1 implies that the pseudosquares prime test takes arithmetic operations to test integers for primality, given a table of pseudosquares If . has no prime divisors below then pseudosquares up to , suffice. See Reference 21Reference 28.
So, under the assumption of Conjecture 1, the cost of applying the pseudosquares prime test to all the integers after they all pass a base-2 pseudoprime test is at most proportional to
and this is for .
The space used is dominated by the length of the sieve intervals and the space needed to store the primes in which is , bits.
This completes the proof of Theorem 2.
As mentioned previously, we implemented several versions of our second algorithm to see what we could compute. We looked for new computational records that were within reach of our university’s nice but aging hardware. Below we discuss some of the results of those computations. Some of the implementation details are specific to a particular computation. Here are four remarks about implementation details that these computations had in common.
We wrote our programs in C++ using the GNU compiler under Linux. GMP was used for multiprecision arithmetic when necessary. Note that it was fairly easy to write the code such that GMP was needed only on occasion and for prime tests.
We used MPI and ran our code on Butler University’s cluster Big Dawg. This machine has 16 compute nodes with 12 cores (2 CPUS), each at optimal capacity. Our average utilization was around 150 of the 192 cores due to compute nodes going down from time to time. The CPU is the Intel Xeon CPU E5-2630 0 @ 2.30GHz with 15 MB cache, with 6 cores per CPU.
To parallelize the algorithm, we striped on the residues In other words, all . processors stepped through all the residues but only sieved every residue for primes. This meant there was very little communication overhead, except for when periodic checkpoints were done, about every 15-30 minutes. th
We usually chose our wheel size ( and sieve intervals so that the size of each interval )( was at most a few megabytes so that it would fit in the CPU cache. We used a )vector<bool>, which packs bits.
For larger values of we observed that when sieving by smaller primes , by each of the residues, we might find that almost all the bits of the current interval were cleared long before we reached the sieving limit so we created a simple early-abort strategy that was able to save time. ,
The very few remaining bits were tested with the base-2 strong pseudoprime test even though we had not sieved all the way to We also, then, replaced the use of the pseudosquares prime test with strong pseudoprime tests .Reference 22 using results from Reference 30 so that only a few bases were needed, due to the spotty trial-division information.
We found that, especially for larger our algorithm spent more time on sieving than prime testing. As mentioned previously, for , the prime testing dominates the running time in practice, and it is worthwhile to use so that prime testing is not required.
4.1. Twin primes and Brun’s constant
Let count the twin prime pairs with and let be the sum of the reciprocals of their elements. Thomas Nicely computed these functions up to (see http://www.trnicely.net/\#PI2X). We verified his computations to 14 digits and extended the results to A portion of our computational results are in the table below. .
10304195697298 1.8304844246583 19831847025792 1.8318080634324 29096690339843 1.8325599218628 38196843833352 1.8330837014776 47177404870103 1.8334845790134 56064358236032 1.8338086822020 64874581322443 1.8340803303554 73619911145552 1.8343139034256 82309090712061 1.8345186031523 90948839353159 1.8347006694414
The last section of Klyve’s PhD thesis Reference 18 describes how to use this information to derive bounds for Brun’s constant.
We have four remarks on our algorithm implementation:
As mentioned above, for small like it is more efficient to set , so that sieving also determines primality, thereby avoiding base-2 strong pseudoprime tests and primality tests.
We computed using Kahan summation Reference 17 with the long double data type in C++, which gave us 17 digits, 14 of which were accurate; Thomas Nicely has data with 53 digits of accuracy. The partial sums were accumulated in 10,000 buckets for each process, and then the buckets were in turn added up across processes using Kahan summation.
Our computation took roughly 3 weeks of wall time, which included at least one restart from a checkpoint. Our verification of Nicely’s work to took 42 hours.
We used a wheel with Note that this is roughly . There were . residues to sieve.
See OEIS.org sequence A007508, the number of twin prime pairs below .
4.2. Quadruplet primes
A related sum involves the reciprocals of the elements of the prime tuple Let . count these tuplets up to and let , be the sum of the reciprocals of their elements. Thomas Nicely computed these functions up to We extended this computation and partial results are in the table below. The first two lines are Thomas Nicely’s own results, which we verified. .
25379433651 0.8704776912340 46998268431 0.8704837109481 67439513530 0.8704870310432 87160212807 0.8704893020026 106365371168 0.8704910169467 125172360474 0.8704923889088 143655957845 0.8704935288452 161868188061 0.8704945017556 179847459283 0.8704953489172 197622677481 0.8704960981105
This computation took about 4 days, and we used a separate program rather than looking for pairs of twin primes in the first program. Even though is large enough to use prime tests, we found that sieving to was faster in practice.
We used a wheel with which gave , residues.
See OEIS.org sequence A050258, the number of prime quadruplets with largest member .
4.3. Cunningham chains
We have two new computational results for Cunningham chains.
We found the smallest chain of length 15 of the first kind, and it begins with the prime
The next four chains of this length of the first kind begin with
This computation took roughly a month of wall time. Here we used wheel size with , residues to sieve. Note that is a rather badly behaved prime for larger Cunningham chains (only residues are excluded), so we left it out of the wheel.
See OEIS.org sequence A005602, smallest prime beginning a Cunningham chain of length (of the first kind).
In 2008 Jaroslaw Wroblewski found a Cunningham chain of length 17 of the first kind, starting with
and we were able to show that this is in fact the smallest such chain of that length.
This computation took roughly three months of wall time. We used with , residues to sieve. With roughly three times as many residues as the previous computation, it took roughly three times as long to complete.
5. Discussion and future work
In summary, we have described and analyzed two algorithms for finding primes in patterns and then shown that the second of these algorithms is quite practical by performing a few computations.
We have some ideas for future work.
In the Introduction, we mentioned that our algorithms could be used to find Carmichael numbers by finding prime triplets that satisfy the pattern , , but we have not yet done that computation ,Reference 6.
A natural extension to our algorithms here is to allow the linear polynomials to potentially be higher degree, irreducible polynomials. See Schinzel’s Hypothesis H. (See Reference 27 and Reference 7, §1.2.2) and the Bateman-Horn conjecture Reference 4.)
An algorithm for twin primes with space roughly that runs in time would be nice.
We wish to thank Frank Levinson for his generous gift that funds the Butler cluster supercomputer Big Dawg and the Holcomb Awards Committee for their financial support of this work.
Thanks to Thomàs Oliveira e Silva for several helpful comments on an earlier version of this paper.
We also wish to thank two anonymous referees for many detailed, helpful comments on earlier versions of this paper. In particular, one of the referees contributed the analysis of the implied constant in §2.3.