Choosing primes uniformly at random

I'm interested in efficient methods of generating prime numbers in a given range [a, b] (or with a given number of bits/digits, etc.). By "efficient" I mean minimizing time, randomness, and perhaps other resources in the average case. For concreteness, let's say the prime is to be selected from the range $(2^n,\,2^{n+1})$.

The naive method is to generate numbers in the range at random and test for primality. This takes $n$ bits per iteration with an average of about $(n+1)\log 2-1\approx n\log 2$ iterations. In total this takes $n^2\log 2$ random bits and $n\log 2$ primality tests.

An obvious improvement is to use only odd numbers (the special case of 2 can be handled separately with O(1) overhead), reducing random bits to $(n-1)\frac{(n+1)\log 2-1}{2}\approx \frac{\log 2}{2}n^2$. This trick can be extended: generate a number in [a/30, b/30], multiply by 30, and choose the residue mod 30 with another 3 bits. (The ends can, once again, be handled with only O(1) extra work.) This reduces the number of bits per iteration by only a small amount (1.9 in this case) but reduces the number of iterations by a factor of $30/\varphi(30)=3.75$ over the naive method (87% improvement over generating just odds).

The extreme version of this algorithm would be to compute $\pi(b)$ and $\pi(a-1)$ and choose an index at random; this takes the information-theoretic minimum number of bits, about $n-\lg n+O(1)$, at the cost of a truly atrocious runtime, something worse than $2^{n/2}$ with current algorithms.

Are there any known algorithms with a better time/randomness tradeoff than this family of algorithms, but which still generates the primes uniformly? (Algorithms like nextprime(random(...)) which are not uniform are excluded, even if, like Joye-Paillier, they are close to uniform.)

• Does this algorithm need to work for any $n\in\Bbb{N}$? If it only needs to work for "reasonably" large (not arbitrarily large) $n$, you could simply select from a list of all primes up to a certain point. It's "cheating" (and potentially memory-heavy; I don't know how many primes are known), but if you're interested in practicality... – KSmarts Dec 15 '14 at 21:38
• I don't know anything about this really. Just one thought/suggestion. Have you checked out how the RSA private key generation algorithms do this? Granted, key generation may not be very time critical there. – Jyrki Lahtonen Dec 15 '14 at 21:38
• @KSmarts: I'm happy to let you pick from a list, but I'm not allowing precomputation (and need arbitrary $n$), so the cost of generating the list would need to be included in the algorithm's time. This makes it less efficient (by a factor of ~$2^{n/2}$) than my 'extreme' algorithm. Practically speaking I'd like algorithms to take 'reasonable' space but I'm willing to consider unreasonable algorithms just to get a different perspective. – Charles Dec 15 '14 at 21:45
• "Close to uniform" is a doubtful claim about nextprime(random()). The chance of picking $p\equiv -1\pmod 6$ in $[10^7,5\cdot 10^7]$ with this is $0.51294$ instead of $0.499996$, which is slightly significant (and a consequence of there not being anyprime $\equiv 3\pmod 6$) – Hagen von Eitzen Dec 16 '14 at 7:15

This just adds some information about things you've already mentioned -- I'm quite interested as well in any methods. I note:

To the authors’ knowledge, the only known prime generation algorithms for which the statistical distance to the uniform distribution can be bounded are the one proposed by Maurer [19,20] on the one hand, and the trivial algorithm (viz. pick a random odd integer in the desired interval, return it if it is prime, and try again otherwise) on the other hand.

For n-digit primes, I call my generic range function with the endpoints. That is, I haven't seen any particularly clever methods for exploiting this. If we want to sacrifice some uniformity perhaps we could do something like Fouque/Tibouchi, generating $n-\delta$ digits then looping just constructing $\delta$ digits. Perhaps exploiting some of the simple divisibility rules to quickly cut out multiples of small primes, but that may be not much help over a decent isprime pretest.

As you mention, the random index selection using prime count method has some nice properties, but don't scale. I found it to be faster than other methods up to 18 bits / 6 digits, but the crossover is very implementation specific. I use $${\rm nth\_prime}( \pi({low}) + {\rm irand}(\pi({low},{high})-1) )$$ where the prime counts can be cached if desired (useful for n-bit and n-digit). For small sizes like this my nth prime routine is using a cached segment of 1-byte-per-30 bits and just counts bits (by word with popcount when possible). Speeding this up with some tables might move the crossover a couple more bits but really doesn't seem worth it. For some practical comparison at larger sizes, using this method for 14-digit random primes takes my program about 1.4 seconds vs. 80 microseconds for the trivial method. Regardless of what we do, this method won't work for 100+ bits (the current record for prime count being about $2^{86}$ using many days on a 600+ node cluster).

For n-bit, the best paper I've seen is Fouque and Tibouchi (ePrint 2011 or ICALP 2014). It has quite a bit of discussion of the issue as well as algorithms. Their algorithms don't set the ${n}$th bit so need some modification to generate in the range ($2^n$,$2^{n+1}$). I use a modification of their A1 that always sets the top bit, and checks divisibility by primes to 19 using just native-int calculations on the 63 bits that are generated per iteration. I did not find A2 to run faster, but it would use fewer random bits.

For generic ranges including n-digit, I use index selection for 18 or fewer bits, the trivial method for 64-bit ranges, and for larger ranges I divide the range into partitions of size < $2^{64}$, randomly select a partition, then do the trivial algorithm within the partition. The astute observer will note this skews toward the lower partitions (which contain more primes). It's quite a small effect for n-bit or n-digit, though much larger for ranges like ($0$,$2^{100}$). Definitely not uniform.

I did have an idea for how to address that issue, with lots of handwaving involved, and is just "closer to uniform". Set ${pclow},{pchigh}$ to the lower and upper prime count limits of low and high respectively. This has a small error at large sizes (Axler 2014 or Dusart 2010) and is fast. Generate a uniform random value between $pclow$ and $pchigh$ and compute an nth prime approximation (e.g. 2nd order Cipolla 1902 with third order correction). Repeat if this is outside the actual low/high range (should be unlikely with a decent size range). Now do a trivial algorithm inside a range around this point. The first portion should be accounting for the overall distribution of primes in the range.