1 files changed, 23 insertions, 33 deletions

diff --git a/ass3/q1/twin_prime.py b/ass3/q1/twin_prime.py
index f31f0fb..54fa742 100644
--- a/ass3/q1/twin_prime.py
+++ b/ass3/q1/twin_prime.py
@@ -7,38 +7,34 @@ def even(n):
     return n % 2 == 0
 
 
-def naive_prime_test(n):
-    if n == 1 or even(n):
-        return False
-    for i in range(n - 1, 2, -1):
-        if n * i == 0:
-            return False
-    return True
-
-
-def is_prime(x, k=20):
+# Use the Miller-Rabin algorithm to test primality of x
+def is_prime(x, k):
+    # Base cases
     if even(x):
         return False
     if x in [2, 3]:
         return True
+    # Determine s and t
     s = 0
     t = x - 1
     while even(t):
         s += 1
         t //= 2
+    # Gather witnesses
     witnesses = [random.randrange(2, x - 1) for _ in range(k)]
     for w in witnesses:
-        if power_mod(w, x - 1, x) != 1:
+        if power_mod(w, x - 1, x) != 1:  # Basic Fermat test
             return False
-        for i in range(1, s + 1):
+        for i in range(1, s + 1):  # Advanced test if Fermat fails
             if power_mod(w, 2 ** i * t, x) == 1 and power_mod(w, 2 ** (i - 1) * t, x) not in [1, x - 1]:
                 return False
-    return True
+    return True  # All witnesses must agree to get here
 
 
+# Use repeated-squaring method to determine: base^power % mod
 def power_mod(base, power, mod):
     binary = bin(power)
-    squares = [base % mod]
+    squares = [base % mod]  # Start with base case
     result = squares[0] if binary[-1] == "1" else 1
     for i in range(1, math.floor(math.log2(power)) + 1):
         squares.append(squares[i - 1] ** 2 % mod)
@@ -47,40 +43,34 @@ def power_mod(base, power, mod):
     return result % mod
 
 
+# Generate a twin prime of bit length m. Only tests values of 6x±1.
 def generate_twin_prime(m):
     twin = False
     while not twin:
         start, end = bit_range(m)
-        n = random.randrange(start + 6 - (start % 6), end, 6)
-        # n = random.randrange(2 ** (m - 1) + 6 - ((2 ** (m - 1)) % 6), 2 ** m, 6)
-        # n += 6 - (n % 6)
-        assert n % 6 == 0
-        iterations = max(4, math.ceil(math.log(n, 4)))
-        twin = is_prime(n - 1, iterations) and is_prime(n + 1, iterations)
-    if not (naive_prime_test(n - 1) and naive_prime_test(n + 1)):
-        raise Exception(f"{n - 1}: {naive_prime_test(n - 1)}, {n + 1}: {naive_prime_test(n + 1)}")
-    assert n + 1 < 2 ** m
-    # print(f"{bin(2**m -1)[2:]}\n{bin(n-1)[2:]}\n{bin(n+1)[2:]}")
+        n = random.randrange(start + 6 - (start % 6), end, 6)  # Only test m bit ints that are multiples of 6
+        num_witnesses = max(4, math.ceil(math.log(n, 4)))  # Based on probability 4^{-k} that a witness will be wrong.
+        twin = is_prime(n - 1, num_witnesses) and is_prime(n + 1, num_witnesses)  # Check twin on either side of 6x
     return n - 1, n + 1
 
+
+# Tuple containing range for bit length m
 def bit_range(m):
     return 2 ** (m - 1), 2 ** m
 
-def main():
-    assert len(sys.argv) >= 2
-    twin = generate_twin_prime(int(sys.argv[1]))
-    write_twin(twin)
-    print(twin)
-
 
+# Write a twin in the correct format to file
 def write_twin(twin):
     with open("output_twin_prime.txt", "w") as file:
         file.write(f"{twin[0]}\n{twin[1]}")
 
 
-for i in range(3, 32):
-    print(i)
-    generate_twin_prime(i)
+def main():
+    assert len(sys.argv) >= 2
+    twin = generate_twin_prime(int(sys.argv[1]))
+    write_twin(twin)
+    print(twin)
+
 
 if __name__ == "__main__":
     main()