## Generating Random Large Primes

### March 31, 2017

Sometimes you need to have a large prime, for testing, cryptography, or some other purpose. I’m talking about primes of several hundred to a few thousand digits that are *certified* — proven to be prime — rather than *probable* primes according to a Miller-Rabin or other probabilistic test. Henry Pocklington’s Criterion, which dates to 1914, gives us a way to find such primes quickly. Paulo Ribenboim explains it thus:

Let

pbe an odd prime, and letkbe a positive integer, such thatpdoes not dividekand 1 <k< 2(p+ 1). LetN= 2kp+ 1. Then the following conditions are equivalent:

Nis a prime.- There exists an integer
a, 1 <a<N, such thata^{(N−1)/2}≡ −1 (modN) and gcd(a+ 1,^{k}N) = 1.

This gives us an algorithm for generating large certified primes. Choose *p* a certified prime. Choose 1 ≤ k < 2 × *p* at random; we ignore the last two possibilities for *k*, so we don’t have to worry about *k* being a multiple of *p*. Compute *N*. For each *a* ∈ {2, 3}, test the conditions for primality. If you don’t find a prime, go back and choose a different random *k*. Once you have a prime *N*, “ratchet up” and restart the process with the new certified prime *N* as the *p* of the next step. Continue until *N* is big enough.

Your task is to write a program that generates random large primes. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

Simplistic Python v3.x

Forgot to mention a possible result:

The trial division routine is completely out of order. Sorry! Here’s the corrected procedure:

Here is my favorite simple primality checker:

It uses a prime wheel, so it should be about twice as fast as yours.

@programmingpraxis: Awesome! Thanks you!

Also, I’vve found another bug in my source. The conditions check should read like this:

The GCD bit still bothers me. If k and N are large and of different size, this will take forever and possibly raise a MemoryError.

Your gcd is fine, and will execute quickly; Knuth guarantees it. Most Python programmers would write it like this:

@Millbrae: By the way, I blindly copied the prime checker from my wheel factorization program, and it does too much work. Here’s a corrected version:

And here is the corresponding function that factors integers, which is what I blindly copied from:

All these functions have appeared on Programming Praxis previously.

Thanks again, PP!

Yet there has to be some faster way to get GCD(a**k + 1, N) checked – if k and N are large. Python and PyPy take forever to finish for number 10**12 < n < 10**18.

[Quote from the description]Henry Pocklington’s Criterion, which dates to 1914, gives us a way to find such primes quickly.[/Quote]

… quickly … Quickly … QUICKLY

@Millbrae: If that gcd isn’t fast enough, you might want to look at this.

@Millbrae: If you do gcd(pow(a, k+1, N), N) there are no memory problems and it runs a lot faster.

Thanks a lot, Paul! I’ve been thinking abuot something like this:

Haven’t tested it though…

Latest code in Python:

Sample output:

I’ve tested N at https://www.alpertron.com.ar/ECM.HTM. Everything’s turned out fine.

Additionally, some numbers. If you manually set p = 4089548364052918515677660474680192167456016670586738083090737083455827, the above code will output…

N is then a 140-digit prime ;)

@Milbrae: It should not take six minutes to generate a random 35-digit prime. You seem to misunderstand the algorithm. Start with a small prime, and use it to generate a larger prime. Then recursively use the larger prime that you just generated to generate an even larger prime. And so on. Continue the “ratcheting” process until the end result is the desired size. My Scheme implementation generates a random 35-digit prime in 15ms, which is four orders of magnitude faster, and generates a 140-digit prime in a little bit less than a second:

Well, then maybe this will suffice…

Sample output:

@Milbrae: Much faster. But the result is 212 digits, not the requested 200. Look at my computation of k-lo and k-hi.

Or slightly changed, showing only the generated prime (and its length)

Sample output: