Safe prime numbers in Python

[Paul Rubin]

Skip Montanaro <skip at pobox.com> writes:
>     def safe_primes(last=1):
>         while True:
>             next = gmpy.next_prime(last)
>             if gmpy.is_prime((next-1)/2):
>                 yield next
>             last = next
> 
> Seems to work fairly well.

This doesn't generate a random prime, since the choice of primes by
next_primes is not uniform.  

A simple, reasonably fast way to generate these primes randomly is

Generate random r, where r==3 (mod 4), i.e. r = 4n+3 for some n
Compute s = (r-1)/2

Do fast crude tests on BOTH r and s for primality, i.e. check for
small prime divisors (3,5,7,...).  If either one has a small divisor
throw both away and generate a new pair.  This lets you discard most
candidate pairs quickly.

Now do a more through test (Miller-Rabin or whatever) to make sure
r,s are really prime.

As mentioned, I have some code around that I'll post when I can find
it.  It's pretty straightforward and can generally find a pair within
a minute or so (maybe less, I don't really remember) on my p3-750.