Generate prime of x decimal digits using bit-oriented prime generator

I've got a question on stackoverflow where somebody asks to generate a random 18-digit prime. Unfortunately, the only prime generator is the one from OpenSSL. This prime generator is however geared towards generating primes for RSA and DSA, which means that the only parameter is the size of the prime in bits.

So the question is: how can I generate a prime of x decimal digits using the bit oriented generator? The requirements are that as much as possible of the range is utilized (i.e. the digit at position x being 1..9) and that the chances of a prime is reasonably well distributed over the range (i.e. identical to generating a prime for a generator that would accept x decimal digits as indicator of the range).

It is of course possible to run the prime generator for any of the bit sizes, and it is possible to ask for a re-run as well.

The range of $18$ digit numbers is from $10^{17}$ to $10^{18}-1$. We have $\log_2(10^{17})\approx 56.473, \log_2(10^{18})\approx 59.795$, so the $18$ digit numbers have from $57$ to $60$ bits. First choose the number of bits in your prime. To have an equal chance of choosing any $18$ digit prime, we need to choose the number of bits with a distribution proportional to the number of primes in the range with that many bits. The function PrimePi(n) counts the number of primes below n.

This leads to the following table:

$$\begin {array} {c|c}\text{n}&\text {PrimePi(n)}\\ \hline\\10^{17}&2.55\cdot10^{15}\\2^{57}&3.65\cdot10^{15}\\2^{58}&7.17\cdot10^{15}\\2^{59}&1.410\cdot 10^{16}\\10^{18}&2.413\cdot 10^{16} \end{array}$$

So you want a $57$ bit prime with probability $\frac {3.65-2.55}{24.13-2.55}\approx 0.051$.

Throw a random number to find how many bits you want. If you are not in one of the end brackets, ask for a random prime of that many bits and return the result. If you are in the $57$ bit bracket, keep generating $57$ bit primes until you get one greater than $10^{17}$. If you are in the $60$ bit bracket, generate them until you get one less than $10^{18}$. This should be rather close to equidistribution across the $18$ digit primes.

• This is absolutely correct, but there are issues! See my answer for another viewpoint. – TonyK Dec 23 '15 at 22:05
• This is what I had in my mind. I didn't know that the PrimePi was so easy to operate with. I'd already calculated the bit ranges, but not yet the PrimePi values. I guess that the $n=60$ will take quite some time though, as it seems to disallow most values. Probably best to go just for $2^{59}$ only in practice. – Maarten Bodewes Dec 23 '15 at 23:40
• $2^{60}\approx 1.15\cdot 10^{18}$, so most of the $60$ bit primes will be in range. $2^{56} \approx 7.2\cdot 10^16$, so about $1/3$ of the $57$ bit primes will be too small. You won't be throwing many away either way. – Ross Millikan Dec 24 '15 at 1:08

Are you using Erik Tews' gensafeprime module? If so, there are two issues:

1. All primes generated by this module have the top two bits set. This makes sense in an RSA context, because we need the product of two $n$-bit primes to be $2n$ bits long, not $2n-1$ bits long (there are more sophisticated ways around this problem, but gensafeprime's method is a good practical solution).
2. This module only generates what it calls safe primes, i.e. primes $p$ such that $\frac12(p-1)$ is prime too. This used to be considered a Good Thing (20+ years ago) for arcane number-theoretical reasons, but I think these days the experts regard it as unnecessary.

So you are not getting randomly distributed primes, whatever you do. Therefore it makes no sense to spend too much effort preserving the non-existent randomness of your source primes. In your place, I would use Ross Millikan's method restricted to $56$-, $57$-, and $58$-bit numbers. If you do this, you will never have to discard out-of-range primes, because $1.5 \times 2^{56}$ is an $18$-bit number.

• I'll try and test the method. I didn't think it was needed but the prime generation I was planning to show is BN_generate_prime_ex. The description of the method doesn't say anything about two bits set - I'll try to test this - and the choice of safe primes is optional (Boolean parameter). – Maarten Bodewes Dec 23 '15 at 22:48
• But then you have BN_is_prime_ex available, so you can just generate random $18$-bit numbers until you find a prime, can't you? – TonyK Dec 24 '15 at 8:56
• Yes, that's basically what Robert said in this answer - so I although I accepted Ross answer (as that's what I asked) I actually used the idea from Robert in my answer to the question on SO. Finally, the answer of Ross showed a nice shortcut. Great site this math site; even more wizards than on the crypto site. – Maarten Bodewes Dec 24 '15 at 10:37

There are $22116397130086627$ primes of $18$ decimal digits. Given a random odd number in the range $[10^{17}, 10^{18}]$, its probability of being prime is about $0.049$. So I wouldn't bother with the OpenSSL prime generator. Just generate pseudo-random integers $x$ with $5\times 10^{16} \le x \le 5 \times 10^{17}-1$ until $2x+1$ is prime: on the average it will take just over 20 tries.

• It's not a direct answer, but I guess I can use BN_rand_range and combine it with the primality test in the same library. Thanks!!! – Maarten Bodewes Dec 24 '15 at 0:30