Safe prime numbers in Python

[Skip Montanaro]

    >> Here's what I came up with using gmpy:
    ...

    Scott> I suspect it might be faster if you looked for the lower prime:

    Scott>      def safe_primes(last=1):
    Scott>          while True:
    Scott>              last = gmpy.next_prime(last)
    Scott>              other = last * 2 + 1
    Scott>              if gmpy.is_prime(other):
    Scott>                  yield other

Yes, you're right.  Starting at 2**1000 I found three safe primes quite
quickly.  Using my version I gave up waiting after seeing one safe prime
float by.

Note that the meaning of the last arg changes when converting to your
formulation.  Suppose I call safe_primes(2**1000).  In my formulation, last
identifies the smallest safe_prime I'm interested in.  In your formulation,
last identifies the smallest prime which defines a larger safe prime.  To
make them the same, your definition would have to change slightly:

    def safe_primes(last=1):
        last = long((last-1)/2)   # assuming future float division...
        while True:
            last = gmpy.next_prime(last)
            other = last * 2 + 1
            if gmpy.is_prime(other):
                yield other