Added binary exponentiaion with respect to modulo (#1428)

* Added binary exponentiaion with respect to modulo

* Added miller rabin: the probabilistic primality test for large numbers

* Removed unused import

* Added test for miller_rabin

* Add test to binary_exp_mod

* Removed test parameter to make Travis CI happy

* unittest.main()  # doctest: +ELLIPSIS   ...

* Update binary_exp_mod.py

* Update binary_exp_mod.py

* Update miller_rabin.py

* from .prime_check import prime_check

Co-authored-by: Christian Clauss <cclauss@me.com>

maths/binary_exp_mod.py

@@ -0,0 +1,28 @@
+def bin_exp_mod(a, n, b):
+    """
+    >>> bin_exp_mod(3, 4, 5)
+    1
+    >>> bin_exp_mod(7, 13, 10)
+    7
+    """
+    # mod b
+    assert not (b == 0), "This cannot accept modulo that is == 0"
+    if n == 0:
+        return 1
+
+    if n % 2 == 1:
+        return (bin_exp_mod(a, n - 1, b) * a) % b
+
+    r = bin_exp_mod(a, n / 2, b)
+    return (r * r) % b
+
+
+if __name__ == "__main__":
+    try:
+        BASE = int(input("Enter Base : ").strip())
+        POWER = int(input("Enter Power : ").strip())
+        MODULO = int(input("Enter Modulo : ").strip())
+    except ValueError:
+        print("Invalid literal for integer")
+
+    print(bin_exp_mod(BASE, POWER, MODULO))

maths/miller_rabin.py

@@ -0,0 +1,50 @@
+import random
+
+from .binary_exp_mod import bin_exp_mod
+
+
+# This is a probabilistic check to test primality, useful for big numbers!
+# if it's a prime, it will return true
+# if it's not a prime, the chance of it returning true is at most 1/4**prec
+def is_prime(n, prec=1000):
+    """
+    >>> from .prime_check import prime_check
+    >>> all(is_prime(i) == prime_check(i) for i in range(1000))
+    True
+    """
+    if n < 2:
+        return False
+
+    if n % 2 == 0:
+        return n == 2
+
+    # this means n is odd
+    d = n - 1
+    exp = 0
+    while d % 2 == 0:
+        d /= 2
+        exp += 1
+
+    # n - 1=d*(2**exp)
+    count = 0
+    while count < prec:
+        a = random.randint(2, n - 1)
+        b = bin_exp_mod(a, d, n)
+        if b != 1:
+            flag = True
+            for i in range(exp):
+                if b == n - 1:
+                    flag = False
+                    break
+                b = b * b
+                b %= n
+            if flag:
+                return False
+            count += 1
+    return True
+
+
+if __name__ == "__main__":
+    n = abs(int(input("Enter bound : ").strip()))
+    print("Here's the list of primes:")
+    print(", ".join(str(i) for i in range(n + 1) if is_prime(i)))