All the files for my blog article "An Intermediate Guide To RSA":

https://dusty.phillips.codes/2018/09/13/an-intermediate-guide-to-rsa/

.gitignore

@@ -0,0 +1,4 @@
+.pytest_cache
+.python-version
+.vscode
+__pycache__

cy_is_prime.pyx

@@ -0,0 +1,36 @@
+'''
+Attempted to make trial_division faster using cython. It got a 10x
+speedup, but of course, it only works with 64 bit integers. When
+working with 1024 bit or higher RSA keys, it can't handle the math.
+
+To build this file on Linux:
+
+    cython --annotate cy_is_prime.pyx
+    gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing \
+      -I/home/dusty/.pyenv/versions/3.7.0/include/python3.7m \
+      -o cy_is_prime.so \
+      cy_is_prime.c
+
+'''
+
+#cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True
+import math
+import random
+
+
+def trial_division(long long py_number):
+    cdef long long number = py_number
+    cdef long long redundant_above = math.ceil(math.sqrt(number))
+    if number <= 1:
+        return 0
+    if number == 2 or number == 3 or number ==5:
+        return 1
+    if not number % 2 or not number % 3 or not number % 5:
+        return 0
+
+    cdef long long i = 5
+    while i < redundant_above +1:
+        if not number % i or not number % (i + 2):
+            return 0
+        i += 6
+    return 1

is_prime.py

@@ -0,0 +1,51 @@
+import math
+import random
+
+FIRST_PRIMES = [2, 3, 5, 7, 13, 17, 19, 23]
+
+
+def trial_division(number: int) -> int:
+    if number <= 1:
+        return False
+    for prime in FIRST_PRIMES:
+        if number == prime:
+            return True
+        if not number % prime:
+            return False
+
+    redundant_above = math.ceil(math.sqrt(number))
+    for i in range(5, redundant_above + 1, 6):
+        if not number % i or not number % (i + 2):
+            return False
+    return True
+
+
+def miller_rabin(number: int) -> int:
+    if number <= 1:
+        return False
+    for prime in FIRST_PRIMES:
+        if number == prime:
+            return True
+        if not number % prime:
+            return False
+
+    odd_factor = number - 1
+    mod_levels = 0
+    while odd_factor % 2 == 0:
+        odd_factor = odd_factor // 2
+        mod_levels += 1
+
+    for trials in range(40):
+        witness = random.randrange(2, number - 1)
+        mod_pow = pow(witness, odd_factor, number)
+        if mod_pow == 1:
+            continue
+        for i in range(mod_levels):
+            if mod_pow == (number - 1):
+                break
+            else:
+                i = i + 1
+                mod_pow = (mod_pow ** 2) % number
+        else:
+            return False
+    return True

notes.md

@@ -0,0 +1,39 @@
+ *  generate two large random primes. well how does one do that?
+ *  primes cannot be calculated, they can only be confirmed
+ *  so generate two random numbers and check if they are prime
+ *  how to check if prime?
+     *  [There are several algorithms](https://en.wikipedia.org/wiki/Primality_test)
+     *  Tested my implementation by copying https://primes.utm.edu/lists/small/10000.txt into a python set and comparing to mine
+
+I tried this but too slow:
+
+    def is_prime(number):
+        if number <= 1:
+            return False
+        for prime in FIRST_PRIMES:
+            if number == prime:
+                return True
+            if not number % prime:
+                return False
+
+        redundant_above = math.ceil(math.sqrt(number))
+        for i in range(FIRST_PRIMES[-1], redundant_above + 1, 6):
+            if not number % i or not number % (i + 2):
+                return False
+        return True
+
+dicked around with cython a bit, but fuggit, I just pip installed sympy
+
+
+ *  The secrets module is cryptographically secure https://docs.python.org/3/library/secrets.html
+ *  since RSA is usually quoted in "x-bit encryption", try to generate an appropriate number of bits. each of the two numbers should have about half the number of bits
+ *  use namedtuple so it can't be changed
+
+ *  What the hell is a totient? Oh... it's just (p-1)*(q-1) why didn't you say so?
+
+ *  algorithm says "select e" such that. Some say it can be hardcoded to 65537 (2**16)
+ *  can use math.gcd, but implemented the euclidean algorithm for better understanding
+     *  The iterative version is hard to understand but honestly it's just
+     *  see https://brilliant.org/wiki/extended-euclidean-algorithm/
+ *  totient should/could use carmeichl instead of euler totient
+ *  Algorithm for modinv seems convoluted, but it's just "We observed that a number x had an inverse mod 26 (i.e., a number y so that xy = 1 mod 26) if and only if gcd(x, 26) = 1." http://www-math.ucdenver.edu/~wcherowi/courses/m5410/exeucalg.html

rsa.py

@@ -0,0 +1,75 @@
+import typing
+from dataclasses import dataclass
+from secrets import randbits
+
+from .is_prime import miller_rabin as is_prime
+
+
+@dataclass
+class RSAKey:
+    modulus: int
+    pub_exponent: int
+    priv_exponent: int
+
+
+def guarantee_bits(value: int, num_bits: int) -> int:
+    """
+    Guarantees top and bottom bit are set 1.
+    set top bit so that it is exactly the right number of bits
+    bottom bit since all primes are odd, so may as well ignore them
+    from the start
+    """
+    return value & (1 << num_bits) - 1 | (1 | 1 << (num_bits - 1))
+
+
+def exact_randbits(num_bits: int) -> int:
+    return guarantee_bits(randbits(num_bits), num_bits)
+
+
+def random_prime(num_bits: int) -> int:
+    number = exact_randbits(num_bits)
+    while not is_prime(number):
+        number = exact_randbits(num_bits)
+    return number
+
+
+def _xgcd(b: int, a: int) -> typing.Tuple[int, int]:
+    "returns (gcd, multiple of modinv)"
+    x0, x1 = 1, 0
+    while a != 0:
+        x0, x1 = x1, x0 - b // a * x1
+        b, a = a, b % a
+    return b, x0
+
+
+def gcd(b: int, a: int) -> int:
+    return _xgcd(b, a)[0]
+
+
+def modinv(base: int, modulus: int) -> int:
+    gcd, x = _xgcd(base, modulus)
+    if gcd != 1:
+        raise ValueError(f"No modular inverse for {base} {modulus}")
+    return x % modulus
+
+
+def generate_key(num_bits: int) -> RSAKey:
+    prime1 = random_prime(num_bits // 2)
+    prime2 = random_prime(num_bits // 2)
+    while prime2 == prime1:
+        prime2 = random_prime(num_bits // 2)
+    key_modulus = prime1 * prime2
+    totient = (prime1 - 1) * (prime2 - 1)
+    pub_exponent = random_prime(num_bits - 2)
+    while pub_exponent >= totient or gcd(pub_exponent, totient) != 1:
+        pub_exponent = random_prime(num_bits - 2)
+    priv_exponent = modinv(pub_exponent, totient)
+    return RSAKey(key_modulus, pub_exponent, priv_exponent)  # noqa: T484
+
+
+def encrypt(key: RSAKey, message: int) -> int:
+    return pow(message, key.pub_exponent, key.modulus)
+
+
+def decrypt(key: RSAKey, ciphertext: int) -> int:
+    return pow(ciphertext, key.priv_exponent, key.modulus)