1 files changed, 98 insertions, 0 deletions
diff --git a/ass3/q1/twin_prime.py b/ass3/q1/twin_prime.py
new file mode 100644
index 0000000..aaeeef8
--- /dev/null
+++ b/ass3/q1/twin_prime.py
@@ -0,0 +1,98 @@
+import random
+import math
+
+
+def even(n):
+    return n % 2 == 0
+
+
+def is_really_prime(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 pi2(n):
+    c = 0.661618
+    return 2 * c * (n / (math.log(n) ** 2))
+
+
+def is_prime(x, k=20):
+    if even(x):
+        return False
+    if x in [2, 3]:
+        return True
+    s = 0
+    t = x - 1
+    while even(t):
+        s += 1
+        t //= 2
+    assert math.floor(t) == t
+    witnesses = [random.randrange(2, x - 1) for _ in range(k)]
+    for w in witnesses:
+        if power_mod(w, x - 1, x) != 1:
+            return False
+        for i in range(1, s + 1):
+            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 all(map(lambda w: is_prime(x, w), [2, 3, 5, 7, 11, 13, 17]))
+
+
+def power_mod(base, power, mod):
+    binary = bin(power)
+    squares = [base % mod]
+    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)
+        if binary[-i - 1] == "1":
+            result *= squares[i]
+    return result % mod
+
+
+# for i in range(3, 100):
+#     print(i, is_prime(i), is_really_prime(i))
+#     if is_prime(i) != is_really_prime(i):
+#         print("Liar", i, is_prime(i), is_really_prime(i))
+
+# is_really_prime(9)
+
+
+# for a in filter(is_really_prime, range(50000, 500000)):
+#     if not miller_rabins(a, max(4, math.ceil(math.log(a, 4)))):
+#         raise Exception(f"Didn't work for {a}")
+
+
+def generate_twin_prime(m):
+    twin = False
+    while not twin:
+        n = random.randrange(2 ** (m - 1), 2 ** m)
+        iterations = max(4, math.ceil(math.log(n, 4)))
+        twin = is_prime(n - 1, iterations) and is_prime(n + 1, iterations)
+    if not (is_really_prime(n - 1) and is_really_prime(n + 1)):
+        raise Exception(f"{n-1}: {is_really_prime(n - 1)}, {n+1}: {is_really_prime(n + 1)}")
+    return n - 1, n + 1
+
+
+for bits in range(3, 64):
+    print(f"==== {bits} bits ==== {2 ** (bits - 1), 2 ** bits - 1}")
+    for _ in range(1000):
+        generate_twin_prime(bits)
+        if _ % 50 == 0:
+            print(".", end="")
+    print()
+
+# i = 0
+# while True:
+#     i += 1
+#     base, power, mod = random.randint(1, 1000), random.randint(1, 5000), random.randint(1, 1000)
+#     mine = power_mod(base, power, mod)
+#     real = (base ** power) % mod
+#     pass_test = mine == real
+#     if i % 10000 == 0:
+#         print(i)
+#     if not pass_test:
+#         # print(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.")
+#         raise Exception(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.")