Random RSA/miller-rabin/euclidean

py/milrab.py

@@ -0,0 +1,78 @@
+import random
+import time
+
+
+def lcg(seed, a=1664525, c=1013904223, m=2**32):
+    while True:
+        seed = (a * seed + c) % m
+        yield seed  # Produces an infinite stream of numbers
+
+
+rand_gen = lcg(42)
+
+
+def randint(a, b) -> int:
+    return next(rand_gen) % (b - a + 1) + a
+
+
+def milrab(n: int, iter: int = 5) -> bool:
+    "Miller-Rabin probabilistic primality test"
+
+    # These should need no further introduction
+    if n < 2:
+        return False
+    if n in (2, 3):
+        return True
+    if n % 2 == 0:
+        return False
+
+    # Keep in mind that here we know that m will be an even number
+    m = n - 1
+    k = 0
+
+    # Executed at least once
+    while m % 2 == 0:
+        m //= 2  # Int divide for int result type
+        k += 1
+
+    for _ in range(iter):
+        a = randint(2, n - 2)  # Any integer in range 1 < a < n - 1
+        b = pow(a, m, n)
+
+        # Valid only for the first iteration
+        if b == 1 or b == n - 1:
+            continue  # Likely a prime
+
+        # The rest, we start from 0
+        for _ in range(k - 1):
+            b = pow(b, 2, n)
+            if b == n - 1:
+                break  # May be prime
+            if b == 1:  # For sure not prime
+                return False
+
+        else:  # Python way to branch if loop never calls break
+            return False
+
+    return True
+
+
+assert milrab(53, 1) == True
+assert milrab(53, 2) == True
+assert milrab(53, 3) == True
+assert milrab(53, 4) == True
+assert milrab(203, 10) == False
+
+assert milrab(10, 4) == False
+
+start = time.perf_counter()
+primes = []
+for p in range(2**15, 2**16):
+    if milrab(p, 10):
+        primes.append(p)
+end = time.perf_counter()
+print(f"Time taken: {end - start:.6f} seconds")
+
+
+for k in range(10):
+    assert milrab(53, k) == True

py/rsa.py

@@ -0,0 +1,95 @@
+from Crypto.Util.number import getPrime, inverse
+import random
+
+
+def miller_rabin(n, k=5):
+    if n == 2 or n == 3:
+        return True
+    if n <= 1 or n % 2 == 0:
+        return False
+
+    # Write n-1 as d * 2^r
+    r = 0
+    d = n - 1
+    while d % 2 == 0:
+        d //= 2
+        r += 1
+
+    # Perform the test k times
+    for _ in range(k):
+        a = random.randint(2, n - 2)
+        x = pow(a, d, n)  # x = a^d % n
+        if x == 1 or x == n - 1:
+            continue
+
+        # Keep squaring x and check for n-1
+        for _ in range(r - 1):
+            x = pow(x, 2, n)
+            if x == n - 1:
+                break
+        else:
+            return False
+
+    return True
+
+
+# Euclidean algorithm
+def extended_gcd(a, b):
+    """Computes the GCD of a and b, along with coefficients x and y such that ax + by = gcd(a, b)."""
+    if b == 0:
+        return a, 1, 0
+
+    g, x1, y1 = extended_gcd(b, a % b)
+
+    x = y1
+    y = x1 - (a // b) * y1
+    return g, x, y
+
+
+def mod_inverse(a, m):
+    g, x, _ = extended_gcd(a, m)
+    if g != 1:  # a and m must be coprime
+        raise ValueError(f"{a} has no modular inverse modulo {m}")
+    return x % m
+
+
+def my_prime(bits: int = 1024) -> int:
+    i = random.randint(1 << (bits - 1), 1 << bits)
+    return i if miller_rabin(i) else my_prime(bits)
+
+
+print(my_prime(8))
+
+
+# Key Generation
+def generate_keys(bits: int = 1024) -> tuple[tuple[int, int], tuple[int, int]]:
+    p = my_prime(bits // 2)
+    q = my_prime(bits // 2)
+    n = p * q
+    phi = (p - 1) * (q - 1)
+    e = 65537  # Common public exponent
+    d = inverse(e, phi)
+    return (e, n), (d, n)  # Public, Private keys
+
+
+# Encryption
+def encrypt(msg: int, pubkey: tuple[int, int]) -> int:
+    e, n = pubkey
+    return pow(msg, e, n)
+
+
+# Decryption
+def decrypt(cipher: int, privkey: tuple[int, int]) -> int:
+    d, n = privkey
+    return pow(cipher, d, n)
+
+
+# Example Usage
+pub, priv = generate_keys()
+pub = (177349751, 2144805071)
+priv = (1698859991, 2144805071)
+message = 42
+ciphertext = encrypt(message, pub)
+plaintext = decrypt(ciphertext, priv)
+
+print(f"Message: {message}, Ciphertext: {ciphertext}, Decrypted: {plaintext}")

py/sketch.py

@@ -0,0 +1,97 @@
+import random
+
+
+def miller_rabin(n, k=5):
+    if n == 2 or n == 3:
+        return True
+    if n <= 1 or n % 2 == 0:
+        return False
+
+    # Write n-1 as d * 2^r
+    r = 0
+    d = n - 1
+    while d % 2 == 0:
+        d //= 2
+        r += 1
+
+    # Perform the test k times
+    for _ in range(k):
+        a = random.randint(2, n - 2)
+        x = pow(a, d, n)  # x = a^d % n
+        if x == 1 or x == n - 1:
+            continue
+
+        # Keep squaring x and check for n-1
+        for _ in range(r - 1):
+            x = pow(x, 2, n)
+            if x == n - 1:
+                break
+        else:
+            return False
+
+    return True
+
+
+assert miller_rabin(0) == False
+assert miller_rabin(1) == False
+assert miller_rabin(2) == True
+assert miller_rabin(3) == True
+assert miller_rabin(4) == False
+
+
+def primitive_isprime(num):
+    for n in range(2, int(num**0.5) + 1):
+        if num % n == 0:
+            return False
+    return True
+
+
+assert not miller_rabin(8)
+assert miller_rabin(11)
+
+
+def mod_inverse(a, m):
+    return pow(a, -1, m)
+
+
+assert mod_inverse(3, 7) == 5
+assert mod_inverse(10, 17) == 12
+assert mod_inverse(7, 13) == 2
+assert mod_inverse(65537, 3445361432) is not None
+
+
+def gcd(x, y):
+    while y:
+        x, y = y, x % y
+    return abs(x)
+
+
+assert gcd(3, 9) == 3
+assert gcd(3, 4) == 1
+
+p = 60737
+q = 56713
+assert miller_rabin(p)
+assert miller_rabin(q)
+
+n = p * q
+phi_n = (p - 1) * (q - 1)
+assert not miller_rabin(n)
+assert not miller_rabin(phi_n)
+
+e = 65537
+assert e == (1 << 16) | 0x1
+assert gcd(e, phi_n) == 1
+
+d = mod_inverse(e, phi_n)
+assert d is not None and (e * d) % phi_n == 1
+
+m = 69
+
+c = pow(69, e, n)
+dec = pow(c, d, n)
+
+print("Ciphertext: ", c)
+print("Decrypted: ", dec)
+
+assert dec == m