diff --git a/py/milrab.py b/py/milrab.py new file mode 100644 index 0000000..be9adab --- /dev/null +++ b/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 diff --git a/py/rsa.py b/py/rsa.py new file mode 100644 index 0000000..5ed0396 --- /dev/null +++ b/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}") diff --git a/py/sketch.py b/py/sketch.py new file mode 100644 index 0000000..025c792 --- /dev/null +++ b/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