# https://ask.sagemath.org/question/44765/randomly-generate-a-safe-prime-of-given-length/
# https://crypto.stackexchange.com/questions/66076/how-to-efficiently-generate-a-random-safe-prime-of-given-length
def rdp(nbits=512):
    while True:
        p = random_prime(2^nbits-1, false, 2^(nbits-1))
        if ZZ((p-1)/2).is_prime():
            return p, (p-1)/2

p, ell = rdp(1024)