count-safe-primes-2048.py

#!/usr/bin/env python

# sudo apt-get install python-gmpy python-crypto
# # or
# easy_install gmpy pycrypto

import gmpy
import math
import Crypto
import Crypto.Random
import Crypto.Util.number as number

def count_safe_primes(bits=2048):
    randfunc = Crypto.Random.new().read
    N = 0
    n = 0
    while True:
        N += 1
        if N % 100 == 0 and n > 0 and N > 0:
            p = float(n) / N
            z = 1.95996
            l = (p + z*z/(2*N) - z * math.sqrt((p*(1-p)+z*z/(4*N))/N))/(1+z*z/N)
            u = (p + z*z/(2*N) + z * math.sqrt((p*(1-p)+z*z/(4*N))/N))/(1+z*z/N)
            c = (u - l) / p
            print "%s %s %.7f %.7f %.7f %.7f"  % (N, n, l, p, u, c)
        q = gmpy.mpz(number.getRandomNBitInteger(bits-1, randfunc))
        q = gmpy.next_prime(q)
        p = 2*q + 1
        if gmpy.is_prime(p):
            n += 1

if __name__ == '__main__':
    count_safe_primes()