from random import randrange, getrandbits def sift_primes(primes_list, pc): for dp in primes_list: if pc % dp == 0: return False primes_list.append(pc) return True def smallprimes(max = 65535): """ Compute a list of known primes. """ primes_list = [2, 3] i = 5 while i <= max: sift_primes(primes_list, i) i += 2 return primes_list def is_prime(n, k=128, def_list=[], min_s=0, totient_coprimes=[65537]): """ Test if a number is prime Args: n -- int -- the number to test k -- int -- the number of combined tests def_list -- int list -- definite known primes, if any, gapless and in ascending order Must start from 2 or 3 min_s -- int -- See below... totient_coprimes -- int list -- What the totient must not be divisible by, which is odd return True if n is prime, and n-1 is divisible by 2^min_s And as long as the probability of non-primeness is extremely small Hence, False could result when n is actually prime This probability can be reduced by setting (min_s = 2), thereby also saving time Setting (min_s > 2) or (min_s = 1) has no known benefits Setting (min_s = -1) will reduce the security, and guarantee that False is only returned, when n is not prime (Not recommended for cryptography) """ # Test if n is not even. # Prevent some error-loops. if min_s < 1 and n == 2: return True if min_s < 2 and n == 3: return True if n == 3: return False if n <= 1 or n % 2 == 0: return False """ First, make sure that the tentative prime - 1 will be coprime with a Public-Key Exponent. """ for dp in totient_coprimes: if dp <= 1 or dp % 2 == 0: continue if (n - 1) % dp == 0: return False # find r and s s = 0 r = n - 1 while r & 1 == 0: s += 1 r //= 2 """ A radical safeguard: Reject if s < 2 Even though n could be prime, Just because the available test would be no stronger than the Fermat Test. """ if min_s > 0 and s < min_s: return False """ Next, check against known primes. """ for dp in def_list: if dp >= n: return True if n % dp == 0: return False if len(def_list) > 0: highest_known_prime = def_list[len(def_list) - 1] highest_known_prime *= highest_known_prime if n < highest_known_prime: return True """ Additional safeguard: If a^r mod n == 1 ... The probability turned out to be 1 / (2^s) . If this probability is exceeded greatly, reject n. """ odd_ones = 0 max_oo = (3 * k) // (1 << (s + 1)) # do k tests for _ in range(k): # Perform Fermat Test first, to save time... a = randrange(2, n - 1) x = pow(a, n-1, n) if x != 1: return False # Now, Proceed with (modified) Miller-Rabin Test. a = randrange(2, n - 1) x = pow(a, r, n) if min_s != -1 and x == 1: odd_ones += 1 if odd_ones > max_oo: return False if x != 1 and x != n - 1: j = 1 while j < s and x != n - 1: x = pow(x, 2, n) if x == 1: return False j += 1 if x != n - 1: return False return True """ Testing above code with non-deterministic assertions, resulting in: n - 1 mod 16 == 0 -> s = 4 n - 1 mod 8 == 0 -> s = 3 But, without using a def_list . If a^r mod n == 1 more than 24 or 48 times, Boom! """ assert is_prime(113, 256) assert is_prime(137, 256) def generate_prime_candidate(length): """ Generate an odd integer randomly Args: length -- int -- the length of the number to generate, in bits return a integer """ # generate random bits p = getrandbits(length) # apply a mask to set MSB and LSB to 1 p |= (1 << length - 1) | 1 return p def generate_prime_number(length=1024): """ Generate a prime Args: length -- int -- length of the prime to generate, in bits return a prime """ p = 4 def_list = smallprimes(65535) # keep generating while the primality test fail while not is_prime(p, 192, def_list, 2): p = generate_prime_candidate(length) return p print(generate_prime_number())