Generating a random prime with a modulus condition in sagemath

Question

I'm wondering if there is a clean way to generate a random prime with a specific modulus condition in sagemath? By modulus condition, I mean, for example, that I might want to generate a random prime that is $1 \pmod{12}$ or $3 \pmod{4}$.

Of course there is random_prime, but I don't see anything in the documentation allowing you to specify a modulus condition. There is a brute force alternative, where you list all the numbers between the desired bounds that satisfy the modulus condition, check if they are prime, and then put all the prime ones in a list and use python functions to choose an element of the list randomly, but I thought perhaps there was a more elegant approach.

Answer 1

It is interesting to look at the code that sage uses for random_prime. It is essentially the following:

def random_prime(lowerbound, upperbound):
    # some input validation that I ignore
    while True:
        p = randint(lowerbound, upperbound)
        if is_prime(p):
            return p

One reason for trying all numbers and then checking primality is that sage wants to give a uniform distribution on all primes in the range. Choosing a random integer and then searching for the smallest prime above it would disproportionately avoid the larger pairs in twin primes, for example --- and sage wants to avoid this.

It is possible to have sage merely look for probable primes, which uses is_pseudoprime in place of is_prime above. Depending on the range, this can be much faster. (And this uses the Baillie-PSW pseudoprimality test, for which there are no currently known exceptions... even though we conjecture that there are infinitely many).

This prompts the potential solution

def random_prime_in_a_mod_b(lowerbound, upperbound, a, b, proof=False):
    """
    Returns a random prime in [lowerbound, upperbound]
    that is congruent to a mod b.
    """
    if not gcd(a, b) == 1:
        raise ValueError("a and b are not coprime")
    if not lowerbound < upperbound:
        raise ValueError("lowerbound is not smaller than upperbound")
    if proof:
        prime_test = is_prime
    else:
        prime_test = is_pseudoprime
    while True:
        n = randint(lowerbound // b, upperbound // b)
        p = n * b + a
        if p < lowerbound or p > upperbound:
            continue
        if prime_test(p):
            return p

This function chooses uniformly among integers congruent to a mod b in a given range. By performing the modulus work ourselves, we don't ever test primes in the wrong congruence class.

Depending on how much you trust the user, it might make sense to use additional input validation, such as checking that a and b are nonnegative, or that the range [lowerbound, upperbound] contains at least one number congruent to a mod b, and so on.