RSA Algorithm: The Foundation of Public-Key Cryptography

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

# Cryptographically Secure Random Number Generation
# ALWAYS use these sources for cryptographic purposes
 
import secrets
import os
 
# Python's secrets module (Python 3.6+)
# This is the recommended approach for cryptographic randomness
 
def generate_random_bits(bit_length: int) -> int:
    """
    Generate a cryptographically secure random integer of specified bit length.
    
    Uses secrets.randbits() which sources entropy from the operating system's
    cryptographic random number generator.
    """
    return secrets.randbits(bit_length)
 
def generate_random_bytes(byte_count: int) -> bytes:
    """
    Generate cryptographically secure random bytes.
    
    Equivalent to reading from /dev/urandom on Unix systems.
    """
    return secrets.token_bytes(byte_count)
 
# Alternative: Direct OS entropy (equivalent to secrets internally)
def generate_random_bytes_os(byte_count: int) -> bytes:
    """
    Read random bytes directly from the OS entropy source.
    
    - Linux/macOS: /dev/urandom
    - Windows: CryptGenRandom
    """
    return os.urandom(byte_count)
 
# For RSA key generation, we need large random odd numbers as prime candidates
def generate_prime_candidate(bit_length: int) -> int:
    """
    Generate a random odd number of exactly the specified bit length.
    
    Requirements:
    1. Must be exactly bit_length bits (MSB and LSB are 1)
    2. Must be odd (LSB = 1) - even numbers > 2 aren't prime
    3. Must come from a CSPRNG
    """
    # Generate random bits
    candidate = secrets.randbits(bit_length)
    
    # Set MSB to ensure exact bit length
    candidate |= (1 << (bit_length - 1))
    
    # Set LSB to ensure odd (even numbers aren't prime)
    candidate |= 1
    
    return candidate
 
# Example usage
if __name__ == "__main__":
    # Generate a 1024-bit prime candidate (half of 2048-bit key)
    candidate = generate_prime_candidate(1024)
    print(f"Prime candidate: {candidate}")
    print(f"Bit length: {candidate.bit_length()}")  # Exactly 1024
    print(f"Is odd: {candidate % 2 == 1}")  # True
 
# WARNING: Never use these for cryptography!
import random  # Uses predictable Mersenne Twister
# random.random(), random.randint(), random.choice() - ALL INSECURE!
 
# The Debian OpenSSL Disaster (2008):
# A developer commented out code that added entropy to OpenSSL's PRNG.
# For 2 years, all SSL keys generated on Debian/Ubuntu systems had only
# 32,767 possible values instead of 2^1024. All were trivially breakable.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

import secrets
 
def miller_rabin_test(n: int, witness: int) -> bool:
    """
    Perform one round of the Miller-Rabin primality test.
    
    Returns True if n passes the test for this witness (probably prime).
    Returns False if n definitely composite.
    
    Mathematical basis:
    For prime p, and p-1 = 2^s * d (d odd), for any witness a:
    Either: a^d ≡ 1 (mod p)
    Or: a^(2^r * d) ≡ -1 (mod p) for some r in [0, s-1]
    """
    if n < 2:
        return False
    if n == 2:
        return True
    if n % 2 == 0:
        return False
    
    # Write n-1 as 2^s * d where d is odd
    s = 0
    d = n - 1
    while d % 2 == 0:
        s += 1
        d //= 2
    
    # Compute a^d mod n
    x = pow(witness, d, n)
    
    # If a^d ≡ 1 (mod n), passes this test
    if x == 1:
        return True
    
    # Check if a^(2^r * d) ≡ -1 (mod n) for any r in [0, s-1]
    for _ in range(s):
        if x == n - 1:  # -1 ≡ n-1 (mod n)
            return True
        x = pow(x, 2, n)
    
    return False  # Definitely composite
 
def is_probably_prime(n: int, iterations: int = 64) -> bool:
    """
    Test if n is probably prime using Miller-Rabin.
    
    With 64 iterations, probability of false positive is < 2^(-128).
    This is astronomically unlikely - more certain than most physical measurements.
    """
    if n < 2:
        return False
    if n == 2 or n == 3:
        return True
    if n % 2 == 0:
        return False
    
    # Test against small primes first (fast rejection of many composites)
    small_primes = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
    for p in small_primes:
        if n == p:
            return True
        if n % p == 0:
            return False
    
    # Perform Miller-Rabin with random witnesses
    for _ in range(iterations):
        witness = secrets.randbelow(n - 3) + 2  # Random in [2, n-2]
        if not miller_rabin_test(n, witness):
            return False
    
    return True
 
def generate_prime(bit_length: int) -> int:
    """
    Generate a random prime number of exactly the specified bit length.
    
    Process:
    1. Generate random odd number with exact bit length
    2. Test primality with Miller-Rabin
    3. Repeat until prime found
    
    Expected iterations: ~0.693 * bit_length / 2 ≈ 355 for 1024 bits
    """
    attempts = 0
    while True:
        attempts += 1
        
        # Generate random odd number with exact bit length
        candidate = secrets.randbits(bit_length)
        candidate |= (1 << (bit_length - 1))  # Set MSB (exact bit length)
        candidate |= 1  # Set LSB (ensure odd)
        
        if is_probably_prime(candidate):
            print(f"Found prime after {attempts} attempts")
            return candidate
 
def generate_rsa_primes(key_bits: int = 2048):
    """
    Generate two suitable primes for RSA key generation.
    
    Requirements for RSA primes:
    1. Each prime is approximately key_bits/2 bits
    2. Primes must be sufficiently different (|p - q| should be large)
    3. Both must be strong primes (resistant to specific attacks)
    """
    prime_bits = key_bits // 2
    
    print(f"Generating first prime ({prime_bits} bits)...")
    p = generate_prime(prime_bits)
    
    print(f"Generating second prime ({prime_bits} bits)...")
    while True:
        q = generate_prime(prime_bits)
        
        # Ensure p ≠ q
        if p == q:
            continue
        
        # Ensure |p - q| is large (primes too close are vulnerable)
        # Rule of thumb: difference should be at least 2^(prime_bits - 100)
        min_diff = 1 << (prime_bits - 100)
        if abs(p - q) < min_diff:
            print("Primes too close, regenerating q...")
            continue
        
        break
    
    return p, q
 
# Example usage
if __name__ == "__main__":
    # Generate 2048-bit RSA primes (this takes several seconds)
    # p, q = generate_rsa_primes(2048)
    
    # For demonstration, use smaller primes
    p = generate_prime(512)
    q = generate_prime(512)
    n = p * q
    
    print(f"p = {p}")
    print(f"q = {q}")
    print(f"n = p * q ({n.bit_length()} bits)")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136

import math
import secrets
 
def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
    """
    Extended Euclidean Algorithm.
    
    Returns (gcd, x, y) where gcd(a,b) = ax + by.
    Used to compute modular multiplicative inverses.
    """
    if a == 0:
        return b, 0, 1
    
    gcd, x1, y1 = extended_gcd(b % a, a)
    x = y1 - (b // a) * x1
    y = x1
    
    return gcd, x, y
 
def mod_inverse(a: int, m: int) -> int:
    """
    Compute the modular multiplicative inverse of a modulo m.
    
    Returns x such that (a * x) % m == 1.
    Requires gcd(a, m) == 1 (a and m must be coprime).
    
    Uses the Extended Euclidean Algorithm.
    """
    gcd, x, _ = extended_gcd(a % m, m)
    
    if gcd != 1:
        raise ValueError(f"Modular inverse doesn't exist: gcd({a}, {m}) = {gcd}")
    
    return (x % m + m) % m  # Ensure positive result
 
def compute_rsa_keys(p: int, q: int, e: int = 65537):
    """
    Compute all RSA key components from primes p and q.
    
    Returns:
        public_key: (e, n)
        private_key: (d, n)
        key_info: dict with all computed values (for educational purposes)
    """
    # Step 1: Compute modulus n = p * q
    n = p * q
    
    # Step 2: Compute Euler's totient φ(n) = (p-1)(q-1)
    phi_n = (p - 1) * (q - 1)
    
    # Step 2b: Also compute Carmichael's totient λ(n) = lcm(p-1, q-1)
    lambda_n = math.lcm(p - 1, q - 1)
    
    # Step 3: Verify e is coprime to φ(n)
    if math.gcd(e, phi_n) != 1:
        raise ValueError(f"e={e} is not coprime to φ(n)={phi_n}")
    
    # Step 4: Compute private exponent d = e^(-1) mod φ(n)
    # Using Euler's totient
    d_euler = mod_inverse(e, phi_n)
    
    # Using Carmichael's totient (produces smaller d)
    d_carmichael = mod_inverse(e, lambda_n)
    
    # Verify: e * d ≡ 1 (mod φ(n))
    assert (e * d_euler) % phi_n == 1
    assert (e * d_carmichael) % lambda_n == 1
    
    # Both d values work for decryption, Carmichael's is smaller/faster
    d = d_carmichael
    
    public_key = (e, n)
    private_key = (d, n)
    
    key_info = {
        'p': p,
        'q': q,
        'n': n,
        'phi_n': phi_n,
        'lambda_n': lambda_n,
        'e': e,
        'd_euler': d_euler,
        'd_carmichael': d_carmichael,
        'n_bits': n.bit_length(),
    }
    
    return public_key, private_key, key_info
 
def validate_key_pair(public_key, private_key):
    """
    Verify that an RSA key pair is valid by testing encryption/decryption.
    """
    e, n = public_key
    d, _ = private_key
    
    # Test with random messages
    for _ in range(10):
        # Random message less than n
        message = secrets.randbelow(n - 1) + 1
        
        # Encrypt: c = m^e mod n
        ciphertext = pow(message, e, n)
        
        # Decrypt: m = c^d mod n
        decrypted = pow(ciphertext, d, n)
        
        assert decrypted == message, "Key pair validation failed!"
    
    return True
 
# Example: Complete RSA key generation
if __name__ == "__main__":
    # For demonstration, use small primes
    # In practice, use 1024-bit primes for 2048-bit RSA
    p = 61
    q = 53
    
    public_key, private_key, info = compute_rsa_keys(p, q)
    
    print("RSA Key Generation Complete")
    print("=" * 50)
    print(f"Primes: p = {info['p']}, q = {info['q']}")
    print(f"Modulus: n = {info['n']} ({info['n_bits']} bits)")
    print(f"Euler's totient: φ(n) = {info['phi_n']}")
    print(f"Carmichael's totient: λ(n) = {info['lambda_n']}")
    print(f"Public exponent: e = {info['e']}")
    print(f"Private exponent (Euler): d = {info['d_euler']}")
    print(f"Private exponent (Carmichael): d = {info['d_carmichael']}")
    print()
    print(f"Public Key: (e={public_key[0]}, n={public_key[1]})")
    print(f"Private Key: (d={private_key[0]}, n={private_key[1]})")
    print()
    
    # Validate
    if validate_key_pair(public_key, private_key):
        print("✓ Key pair validated successfully!")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

def compute_crt_components(p: int, q: int, d: int):
    """
    Compute CRT optimization components for RSA decryption.
    
    These values enable 4x faster decryption using the Chinese Remainder Theorem.
    """
    # CRT exponents
    dP = d % (p - 1)  # d mod (p-1)
    dQ = d % (q - 1)  # d mod (q-1)
    
    # CRT coefficient
    qInv = mod_inverse(q, p)  # q^(-1) mod p
    
    return dP, dQ, qInv
 
def rsa_decrypt_crt(ciphertext: int, p: int, q: int, dP: int, dQ: int, qInv: int) -> int:
    """
    RSA decryption using Chinese Remainder Theorem optimization.
    
    Approximately 4x faster than standard decryption for large keys.
    
    Process:
    1. Compute m_p = c^dP mod p
    2. Compute m_q = c^dQ mod q  
    3. Combine: m = m_q + q * ((qInv * (m_p - m_q)) mod p)
    """
    # Step 1: Compute partial results
    m_p = pow(ciphertext, dP, p)
    m_q = pow(ciphertext, dQ, q)
    
    # Step 2: Combine using CRT
    # h = qInv * (m_p - m_q) mod p
    h = (qInv * (m_p - m_q)) % p
    
    # m = m_q + h * q
    message = m_q + h * q
    
    return message
 
def rsa_decrypt_standard(ciphertext: int, d: int, n: int) -> int:
    """
    Standard RSA decryption (without CRT).
    
    m = c^d mod n
    """
    return pow(ciphertext, d, n)
 
# Performance comparison
import time
 
def benchmark_decryption(iterations=100):
    """
    Compare CRT vs standard RSA decryption performance.
    """
    # Generate 2048-bit key for realistic benchmark
    # Using smaller key for quick demonstration
    p = generate_prime(512)
    q = generate_prime(512)
    n = p * q
    e = 65537
    phi_n = (p - 1) * (q - 1)
    d = mod_inverse(e, phi_n)
    
    # CRT components
    dP, dQ, qInv = compute_crt_components(p, q, d)
    
    # Generate test ciphertexts
    ciphertexts = [pow(secrets.randbelow(n-1) + 1, e, n) for _ in range(iterations)]
    
    # Benchmark standard decryption
    start = time.time()
    for c in ciphertexts:
        rsa_decrypt_standard(c, d, n)
    standard_time = time.time() - start
    
    # Benchmark CRT decryption
    start = time.time()
    for c in ciphertexts:
        rsa_decrypt_crt(c, p, q, dP, dQ, qInv)
    crt_time = time.time() - start
    
    print(f"Standard decryption: {standard_time:.3f}s")
    print(f"CRT decryption: {crt_time:.3f}s")
    print(f"Speedup: {standard_time/crt_time:.2f}x")
 
# SECURITY NOTE: CRT implementation must check for faults!
# Bellcore attack: If a single-bit error occurs during CRT computation,
# the faulty result can be used to factor n.
# Mitigation: Always verify c^e = m mod n after CRT decryption.