# Re: Safe prime numbers in Python

Discussion in 'Python' started by Skip Montanaro, Jan 14, 2004.

1. ### Skip MontanaroGuest

Tim> Does any know of or have Python code (or C or Fortran code wrapped
Tim> as a Python module) for efficiently finding safe prime numbers -
Tim> that is a number p such that p and (p-1)/2 are both prime?

Here's what I came up with using gmpy:

import gmpy

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.

Skip

Skip Montanaro, Jan 14, 2004

2. ### Scott David DanielsGuest

Skip Montanaro wrote:
> Tim> Does any know of or have Python code (or C or Fortran code wrapped
> Tim> as a Python module) for efficiently finding safe prime numbers -
> Tim> that is a number p such that p and (p-1)/2 are both prime?
>
> Here's what I came up with using gmpy:
>
> import gmpy
>
> def safe_primes(last=1):
> while True:
> next = gmpy.next_prime(last)
> if gmpy.is_prime((next-1)/2):
> yield next
> last = next
>

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

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

-Scott David Daniels

Scott David Daniels, Jan 14, 2004

3. ### Skip MontanaroGuest

>> 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

Skip Montanaro, Jan 14, 2004
4. ### Skip MontanaroGuest

Ack... I'm an idiot:

Skip> def safe_primes(last=1):
Skip> last = long((last-1)/2) # assuming future float division...

should be last = (last-1)//2 # doh!

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

Skip

Skip Montanaro, Jan 14, 2004
5. ### Tim ChurchesGuest

On Thu, 2004-01-15 at 03:35, Skip Montanaro wrote:
> >> Here's what I came up with using gmpy:

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

Many thanks to Skip Montanaro and Scott Daniels for these suggestions. I
did not suspect that the gmpy library would have functions for prime
finding or testing. With a few days of processing I'll have more safe
primes than I'll ever need. Next step is to find all quadratic residuals
modulo p (where p is a safe prime). It notice that GMP has a function to
calculate the Legendre symbol (which hopefully gmpy wraps), so it should
be simple. Fab!

--

Tim C

PGP/GnuPG Key 1024D/EAF993D0 available from keyservers everywhere
or at http://members.optushome.com.au/tchur/pubkey.asc
Key fingerprint = 8C22 BF76 33BA B3B5 1D5B EB37 7891 46A9 EAF9 93D0

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.0.7 (GNU/Linux)

iD8DBQBABbpJeJFGqer5k9ARAsTRAJ9u7bojn2YbPH/VhalMwYKnMWkcJQCeJRRh
w6sn3S7SbsZnu+AtQ3UFbqA=
=JxVa
-----END PGP SIGNATURE-----

Tim Churches, Jan 14, 2004
6. ### Paul RubinGuest

Skip Montanaro <> 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.

Paul Rubin, Jan 15, 2004