For what it's worth, here is an implementation of the Miller-Rabin primality test that first tries the smallest 13 primes and then tries some random bases:

For what it's worth, here is an implementation of the Miller-Rabin primality test that first tries the smallest 13 primes and then tries some random bases: def is_probable_prime(n, num_random_trials=5): """ An implementation of the Miller-Rabin primality test. If this function returns False, then n is definitely not prime. If this function returns True, then n is probably prime: the probability of n being composite is less than 1/4^(num_random_trials), and nonzero only if n >= 3317044064679887385961981 > 3.3*10^24 > 1.3*2^81 (by https://arxiv.org/abs/1509.00864). """ assert n == abs(int(n)), 'Need a nonnegative integer, got %s' % n if n < 2: return False if n % 2 == 0: return n == 2 (s, d) = (0, n - 1) while d % 2 == 0: (s, d) = (s + 1, d // 2) assert n - 1 == (2**s) * d and d % 2 == 1 small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41] for i in range(len(small_primes) + num_random_trials): a = small_primes[i] if i < len(small_primes) else random.randrange(small_primes[-1], n - 1) if a >= n: break if not (pow(a, d, n) == 1 or any(pow(a, 2**j * d, n) == n - 1 for j in range(s))): return False return True