Sync

Task/Miller-Rabin-primality-test/Python/miller-rabin-primality-test-1.py

@@ -0,0 +1,78 @@
+import random
+
+_mrpt_num_trials = 5 # number of bases to test
+
+def is_probable_prime(n):
+    """
+    Miller-Rabin primality test.
+
+    A return value of False means n is certainly not prime. A return value of
+    True means n is very likely a prime.
+
+    >>> is_probable_prime(1)
+    Traceback (most recent call last):
+        ...
+    AssertionError
+    >>> is_probable_prime(2)
+    True
+    >>> is_probable_prime(3)
+    True
+    >>> is_probable_prime(4)
+    False
+    >>> is_probable_prime(5)
+    True
+    >>> is_probable_prime(123456789)
+    False
+
+    >>> primes_under_1000 = [i for i in range(2, 1000) if is_probable_prime(i)]
+    >>> len(primes_under_1000)
+    168
+    >>> primes_under_1000[-10:]
+    [937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
+
+    >>> is_probable_prime(6438080068035544392301298549614926991513861075340134\
+3291807343952413826484237063006136971539473913409092293733259038472039\
+7133335969549256322620979036686633213903952966175107096769180017646161\
+851573147596390153)
+    True
+
+    >>> is_probable_prime(7438080068035544392301298549614926991513861075340134\
+3291807343952413826484237063006136971539473913409092293733259038472039\
+7133335969549256322620979036686633213903952966175107096769180017646161\
+851573147596390153)
+    False
+    """
+    assert n >= 2
+    # special case 2
+    if n == 2:
+        return True
+    # ensure n is odd
+    if n % 2 == 0:
+        return False
+    # write n-1 as 2**s * d
+    # repeatedly try to divide n-1 by 2
+    s = 0
+    d = n-1
+    while True:
+        quotient, remainder = divmod(d, 2)
+        if remainder == 1:
+            break
+        s += 1
+        d = quotient
+    assert(2**s * d == n-1)
+
+    # test the base a to see whether it is a witness for the compositeness of n
+    def try_composite(a):
+        if pow(a, d, n) == 1:
+            return False
+        for i in range(s):
+            if pow(a, 2**i * d, n) == n-1:
+                return False
+        return True # n is definitely composite
+
+    for i in range(_mrpt_num_trials):
+        a = random.randrange(2, n)
+        if try_composite(a):
+            return False
+
+    return True # no base tested showed n as composite

Task/Miller-Rabin-primality-test/Python/miller-rabin-primality-test-2.py

@@ -0,0 +1,37 @@
+def _try_composite(a, d, n, s):
+    if pow(a, d, n) == 1:
+        return False
+    for i in range(s):
+        if pow(a, 2**i * d, n) == n-1:
+            return False
+    return True # n  is definitely composite
+
+def is_prime(n, _precision_for_huge_n=16):
+    if n in _known_primes or n in (0, 1):
+        return True
+    if any((n % p) == 0 for p in _known_primes):
+        return False
+    d, s = n - 1, 0
+    while not d % 2:
+        d, s = d >> 1, s + 1
+    # Returns exact according to http://primes.utm.edu/prove/prove2_3.html
+    if n < 1373653:
+        return not any(_try_composite(a, d, n, s) for a in (2, 3))
+    if n < 25326001:
+        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5))
+    if n < 118670087467:
+        if n == 3215031751:
+            return False
+        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7))
+    if n < 2152302898747:
+        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11))
+    if n < 3474749660383:
+        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13))
+    if n < 341550071728321:
+        return not any(_try_composite(a, d, n, s) for a in (2, 3, 5, 7, 11, 13, 17))
+    # otherwise
+    return not any(_try_composite(a, d, n, s)
+                   for a in _known_primes[:_precision_for_huge_n])
+
+_known_primes = [2, 3]
+_known_primes += [x for x in range(5, 1000, 2) if is_prime(x)]