The Cryptographic Underground: How Codebreakers Shaped the Digital Revolution

# Simplified Enigma rotor simulation showing the self-encryption prohibition
class EnigmaRotor:
    def __init__(self, wiring, position=0):
        self.wiring = wiring
        self.position = position
    
    def encode(self, letter_index):
        # Forward pass through rotor
        adjusted_index = (letter_index + self.position) % 26
        encoded = self.wiring[adjusted_index]
        # Enigma constraint: no letter maps to itself
        assert encoded != letter_index, "Enigma never encrypts letter to itself"
        return (encoded - self.position) % 26
    
    def step(self):
        self.position = (self.position + 1) % 26
        return self.position == 0  # True if rotor completed full rotation

# The statistical attack that broke Enigma
def analyze_enigma_weakness(ciphertext, known_plaintext):
    """Demonstrates the contradiction method used at Bletchley Park"""
    contradictions = 0
    for i in range(len(ciphertext)):
        if ciphertext[i] == known_plaintext[i]:
            contradictions += 1
    # If contradictions > 0, this key setting is impossible
    return contradictions == 0

# Simplified simulation of Colossus's statistical attack on Lorenz
import numpy as np
from scipy.stats import chi2

def colossus_chi_test(ciphertext, test_key_stream):
    """Simulate Colossus's chi-squared statistical attack"""
    # XOR ciphertext with test key to get potential plaintext
    potential_plaintext = [c ^ k for c, k in zip(ciphertext, test_key_stream)]
    
    # Count character frequencies
    freq_observed = np.zeros(26)
    for char in potential_plaintext:
        if 0 <= char < 26:  # Valid letter
            freq_observed[char] += 1
    
    # Expected frequencies for German text
    freq_expected = np.array([
        8.17, 1.49, 2.78, 4.25, 17.4, 1.66, 3.01, 4.76, 7.55, 0.27,
        1.21, 3.44, 2.53, 9.78, 2.51, 0.79, 0.02, 7.0, 7.27, 6.15,
        4.35, 0.67, 1.89, 0.03, 0.04, 1.13
    ]) * len(potential_plaintext) / 100
    
    # Calculate chi-squared statistic
    chi_squared = np.sum((freq_observed - freq_expected)**2 / freq_expected)
    
    # Lower chi-squared indicates more "German-like" text
    return chi_squared

# Colossus would test thousands of wheel positions using this method
def colossus_wheel_search(ciphertext, max_positions=1000):
    """Simulate Colossus's systematic wheel position search"""
    best_score = float('inf')
    best_position = 0
    
    for position in range(max_positions):
        test_key = generate_lorenz_keystream(position, len(ciphertext))
        score = colossus_chi_test(ciphertext, test_key)
        
        if score < best_score:
            best_score = score
            best_position = position
    
    return best_position, best_score

# Simplified DES S-box demonstration showing NSA modifications
# These substitution boxes were the heart of the NSA controversy

# Original IBM Lucifer S-box (simplified)
lucifer_sbox = [
    [14, 4, 13, 1, 2, 15, 11, 8, 3, 10, 6, 12, 5, 9, 0, 7],
    [0, 15, 7, 4, 14, 2, 13, 1, 10, 6, 12, 11, 9, 5, 3, 8],
    # ... more rows
]

# DES S-box (NSA modified) - S1 example
des_s1 = [
    [14, 4, 13, 1, 2, 15, 11, 8, 3, 10, 6, 12, 5, 9, 0, 7],
    [0, 15, 7, 4, 14, 2, 13, 1, 10, 6, 12, 11, 9, 5, 3, 8],
    [4, 1, 14, 8, 13, 6, 2, 11, 15, 12, 9, 7, 3, 10, 5, 0],
    [15, 12, 8, 2, 4, 9, 1, 7, 5, 11, 3, 14, 10, 0, 6, 13]
]

def sbox_substitution(input_bits, sbox):
    """Perform S-box substitution (core of DES security)"""
    # Extract row (bits 0,5) and column (bits 1-4)
    row = (input_bits & 0x20) >> 4 | (input_bits & 0x01)
    col = (input_bits & 0x1E) >> 1
    
    return sbox[row][col]

# The key insight: NSA modifications made DES resistant to differential cryptanalysis
def differential_characteristic_test(sbox, input_diff):
    """Test S-box resistance to differential cryptanalysis"""
    output_diffs = {}
    
    for input_val in range(64):  # All possible 6-bit inputs
        output1 = sbox_substitution(input_val, sbox)
        output2 = sbox_substitution(input_val ^ input_diff, sbox)
        output_diff = output1 ^ output2
        
        output_diffs[output_diff] = output_diffs.get(output_diff, 0) + 1
    
    # Lower maximum probability = better resistance
    max_prob = max(output_diffs.values()) / 64.0
    return max_prob

# DES S-boxes have maximum differential probability of 1/4
# This was intentionally designed to resist differential attacks

# RSA implementation showing the mathematical elegance
import random
from math import gcd

def is_prime(n, k=5):
    """Miller-Rabin probabilistic primality test"""
    if n < 2: return False
    if n in [2, 3]: return True
    if n % 2 == 0: return False
    
    # Write n-1 as d * 2^r
    d = n - 1
    r = 0
    while d % 2 == 0:
        d //= 2
        r += 1
    
    for _ in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, d, n)
        if x in [1, n - 1]:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True

def generate_large_prime(bits):
    """Generate a large prime for RSA"""
    while True:
        candidate = random.getrandbits(bits)
        candidate |= (1 << bits - 1) | 1  # Ensure it's odd and has correct bit length
        if is_prime(candidate):
            return candidate

def extended_gcd(a, b):
    """Extended Euclidean Algorithm"""
    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(e, phi):
    """Find modular inverse of e modulo phi"""
    gcd, x, _ = extended_gcd(e, phi)
    if gcd != 1:
        raise ValueError("Modular inverse does not exist")
    return x % phi

class RSAKeyPair:
    def __init__(self, bits=2048):
        # Generate two large primes
        self.p = generate_large_prime(bits // 2)
        self.q = generate_large_prime(bits // 2)
        
        # Calculate n and phi(n)
        self.n = self.p * self.q
        self.phi = (self.p - 1) * (self.q - 1)
        
        # Choose e (commonly 65537)
        self.e = 65537
        if gcd(self.e, self.phi) != 1:
            raise ValueError("e and phi(n) are not coprime")
        
        # Calculate private exponent d
        self.d = mod_inverse(self.e, self.phi)
    
    def encrypt(self, message):
        """RSA encryption: c = m^e mod n"""
        if message >= self.n:
            raise ValueError("Message too large for key size")
        return pow(message, self.e, self.n)
    
    def decrypt(self, ciphertext):
        """RSA decryption: m = c^d mod n"""
        return pow(ciphertext, self.d, self.n)
    
    def sign(self, message):
        """RSA signature: s = m^d mod n"""
        return pow(message, self.d, self.n)
    
    def verify(self, signature, message):
        """RSA signature verification: m =? s^e mod n"""
        return pow(signature, self.e, self.n) == message

# Demonstrate RSA's revolutionary capability
rsa = RSAKeyPair(1024)  # Small key for demo
message = 123456789
encrypted = rsa.encrypt(message)
decrypted = rsa.decrypt(encrypted)
print(f"Original: {message}, Encrypted: {encrypted}, Decrypted: {decrypted}")
print(f"Success: {message == decrypted}")

# Simplified demonstration of the Dual_EC_DRBG backdoor mechanism
# Note: This is educational - do not use for actual cryptography!

class EllipticCurve:
    def __init__(self, p, a, b):
        self.p = p  # Prime modulus
        self.a = a  # Curve parameter a
        self.b = b  # Curve parameter b
    
    def point_add(self, P, Q):
        """Elliptic curve point addition"""
        if P is None: return Q
        if Q is None: return P
        
        x1, y1 = P
        x2, y2 = Q
        
        if x1 == x2:
            if y1 == y2:
                # Point doubling
                s = (3 * x1 * x1 + self.a) * pow(2 * y1, -1, self.p) % self.p
            else:
                return None  # Point at infinity
        else:
            # Point addition
            s = (y2 - y1) * pow(x2 - x1, -1, self.p) % self.p
        
        x3 = (s * s - x1 - x2) % self.p
        y3 = (s * (x1 - x3) - y1) % self.p
        
        return (x3, y3)
    
    def scalar_mult(self, k, point):
        """Scalar multiplication of point by integer k"""
        if k == 0: return None
        if k == 1: return point
        
        result = None
        addend = point
        
        while k:
            if k & 1:
                result = self.point_add(result, addend)
            addend = self.point_add(addend, addend)
            k >>= 1
        
        return result

class DualECDRBG:
    """Simplified Dual_EC_DRBG showing the backdoor mechanism"""
    
    def __init__(self, curve, P, Q, seed):
        self.curve = curve
        self.P = P  # Generator point
        self.Q = Q  # Second point (potentially backdoored)
        self.state = seed
    
    def next_bits(self, num_bits):
        """Generate next random bits - this is where the backdoor operates"""
        # s_{i+1} = x-coordinate of s_i * P
        point = self.curve.scalar_mult(self.state, self.P)
        self.state = point[0] if point else 0
        
        # r_i = x-coordinate of s_i * Q  
        output_point = self.curve.scalar_mult(self.state, self.Q)
        output = output_point[0] if output_point else 0
        
        # Return specified number of bits
        return output & ((1 << num_bits) - 1)

# The backdoor: if Q = d * P for known d, then:
# If you know d and observe output r_i = x(s_i * Q) = x(s_i * d * P)
# You can compute s_i and predict all future outputs

def demonstrate_backdoor():
    """Show how the NSA backdoor in Dual_EC_DRBG works"""
    # Simplified curve parameters (toy example)
    p = 2**255 - 19  # Small prime for demo
    curve = EllipticCurve(p, -1, 1)
    
    # Public parameters
    P = (15112221349535400772501151409588531511454012693041857206046113283949847762202, 
         46316835694926478169428394003475163141307993866256225615783033603165251855960)
    
    # Backdoored Q: Q = d * P where NSA knows d
    d = 12345  # NSA's secret backdoor key
    Q = curve.scalar_mult(d, P)
    
    # Initialize PRNG
    rng = DualECDRBG(curve, P, Q, 98765)
    
    # Generate some "random" bits
    output1 = rng.next_bits(128)
    output2 = rng.next_bits(128)
    
    print(f"Dual_EC_DRBG Output 1: {output1:032x}")
    print(f"Dual_EC_DRBG Output 2: {output2:032x}")
    
    # NSA can predict future outputs using backdoor key d
    print("NSA backdoor allows prediction of all future values")
    
    return output1, output2

# The real controversy: NIST standardized this despite academic warnings
demonstrate_backdoor()

# Demonstration of Bitcoin's core cryptographic components
import hashlib
import struct
from datetime import datetime

class MerkleTree:
    """Bitcoin's Merkle tree implementation for transaction verification"""
    
    def __init__(self, transactions):
        self.transactions = transactions
        self.tree = self.build_tree([self.hash_transaction(tx) for tx in transactions])
    
    def hash_transaction(self, transaction):
        """Double SHA-256 hash (Bitcoin's standard)"""
        return hashlib.sha256(hashlib.sha256(transaction.encode()).digest()).digest()
    
    def build_tree(self, hashes):
        """Build Merkle tree bottom-up"""
        if len(hashes) == 1:
            return hashes[0]
        
        # Pad with duplicate if odd number of elements
        if len(hashes) % 2 == 1:
            hashes.append(hashes[-1])
        
        parent_level = []
        for i in range(0, len(hashes), 2):
            combined = hashes[i] + hashes[i + 1]
            parent_hash = hashlib.sha256(hashlib.sha256(combined).digest()).digest()
            parent_level.append(parent_hash)
        
        return self.build_tree(parent_level)
    
    def get_merkle_proof(self, tx_index):
        """Generate proof that transaction is in the block"""
        hashes = [self.hash_transaction(tx) for tx in self.transactions]
        proof = []
        index = tx_index
        
        while len(hashes) > 1:
            if len(hashes) % 2 == 1:
                hashes.append(hashes[-1])
            
            if index % 2 == 0:  # Left child
                proof.append(('right', hashes[index + 1]))
            else:  # Right child
                proof.append(('left', hashes[index - 1]))
            
            # Move to parent level
            parent_hashes = []
            for i in range(0, len(hashes), 2):
                combined = hashes[i] + hashes[i + 1]
                parent_hash = hashlib.sha256(hashlib.sha256(combined).digest()).digest()
                parent_hashes.append(parent_hash)
            
            hashes = parent_hashes
            index //= 2
        
        return proof

class BitcoinBlock:
    """Simplified Bitcoin block structure"""
    
    def __init__(self, previous_hash, transactions, difficulty=4):
        self.previous_hash = previous_hash
        self.transactions = transactions
        self.merkle_tree = MerkleTree(transactions)
        self.merkle_root = self.merkle_tree.tree
        self.timestamp = int(datetime.now().timestamp())
        self.difficulty = difficulty
        self.nonce = 0
        self.hash = self.mine_block()
    
    def calculate_hash(self):
        """Calculate block hash using proof-of-work"""
        block_data = (
            self.previous_hash +
            self.merkle_root +
            struct.pack('>Q', self.timestamp) +
            struct.pack('>L', self.nonce)
        )
        return hashlib.sha256(hashlib.sha256(block_data).digest()).digest()
    
    def mine_block(self):
        """Proof-of-work mining: find nonce that produces hash with leading zeros"""
        target = b'\x00' * self.difficulty  # Simplified difficulty target
        
        while True:
            block_hash = self.calculate_hash()
            if block_hash[:self.difficulty] == target:
                print(f"Block mined! Nonce: {self.nonce}, Hash: {block_hash.hex()[:16]}...")
                return block_hash
            self.nonce += 1
    
    def verify_transaction(self, tx_index, merkle_proof):
        """Verify a transaction is in this block using Merkle proof"""
        tx_hash = self.merkle_tree.hash_transaction(self.transactions[tx_index])
        current_hash = tx_hash
        
        for direction, sibling_hash in merkle_proof:
            if direction == 'left':
                combined = sibling_hash + current_hash
            else:
                combined = current_hash + sibling_hash
            current_hash = hashlib.sha256(hashlib.sha256(combined).digest()).digest()
        
        return current_hash == self.merkle_root

# Demonstrate Bitcoin's cryptographic verification
transactions = [
    "Alice pays Bob 1 BTC",
    "Carol pays Dave 2.5 BTC", 
    "Eve pays Frank 0.1 BTC",
    "Grace pays Henry 3 BTC"
]

# Create and mine a block
block = BitcoinBlock(b'\x00' * 32, transactions, difficulty=2)

# Verify a transaction using Merkle proof
tx_index = 1  # Carol pays Dave
proof = block.merkle_tree.get_merkle_proof(tx_index)
is_valid = block.verify_transaction(tx_index, proof)
print(f"Transaction verification: {is_valid}")
print(f"Merkle proof size: {len(proof)} hashes (vs {len(transactions)} full block)")

# Demonstration of quantum period finding (heart of Shor's algorithm)
# Note: This is a classical simulation for educational purposes

import numpy as np
from math import gcd
import random

def classical_period_finding(a, N):
    """Classical approach to period finding (exponentially slow)"""
    for r in range(1, N):
        if pow(a, r, N) == 1:
            return r
    return None

def quantum_period_finding_simulation(a, N):
    """Simplified simulation of quantum period finding"""
    # In a real quantum computer, this would use quantum Fourier transform
    # We simulate the measurement outcome
    
    # Find the actual period
    actual_period = classical_period_finding(a, N)
    if actual_period is None:
        return None
    
    # Quantum algorithm would find period with high probability
    # Simulate measurement with some noise
    noise = random.uniform(0.9, 1.1)
    measured_period = int(actual_period * noise)
    
    return measured_period

def shors_algorithm_simulation(N):
    """Simplified simulation of Shor's factorization algorithm"""
    if N % 2 == 0:
        return 2, N // 2
    
    # Step 1: Choose random a < N
    while True:
        a = random.randint(2, N - 1)
        
        # Check if gcd(a, N) > 1 (lucky guess)
        g = gcd(a, N)
        if g > 1:
            return g, N // g
        
        # Step 2: Find period of a^x mod N using quantum algorithm
        period = quantum_period_finding_simulation(a, N)
        
        if period is None or period % 2 != 0:
            continue  # Try again with different a
        
        # Step 3: Use period to find factors
        half_period = period // 2
        factor1 = gcd(pow(a, half_period) - 1, N)
        factor2 = gcd(pow(a, half_period) + 1, N)
        
        if factor1 > 1 and factor1 < N:
            return factor1, N // factor1
        if factor2 > 1 and factor2 < N:
            return factor2, N // factor2
        
        # If this didn't work, try again
        continue

# Demonstrate the quantum threat to RSA
def demonstrate_quantum_threat():
    """Show how quantum computing threatens current cryptography"""
    # Small RSA modulus for demonstration
    p, q = 61, 53
    N = p * q
    
    print(f"RSA modulus N = {p} × {q} = {N}")
    
    # Classical factoring is hard
    print("\nClassical factoring attempt:")
    for i in range(2, int(N**0.5) + 1):
        if N % i == 0:
            print(f"Found factors by trial division: {i} and {N // i}")
            break
    
    # Quantum factoring would be efficient
    print("\nSimulated quantum factoring (Shor's algorithm):")
    factor1, factor2 = shors_algorithm_simulation(N)
    print(f"Quantum computer found factors: {factor1} and {factor2}")
    
    return factor1, factor2

# Post-quantum cryptography example: Lattice-based cryptography
class LWECryptosystem:
    """Learning With Errors - a post-quantum cryptographic primitive"""
    
    def __init__(self, n=10, q=97, sigma=1.0):
        self.n = n  # Dimension
        self.q = q  # Modulus
        self.sigma = sigma  # Error standard deviation
        
        # Generate secret key
        self.s = np.random.randint(0, q, n)
        
        # Generate public key (A, b = As + e)
        self.A = np.random.randint(0, q, (n, n))
        e = np.random.normal(0, sigma, n).astype(int) % q
        self.b = (np.dot(self.A, self.s) + e) % q
    
    def encrypt(self, message_bit):
        """Encrypt a single bit using LWE"""
        # Choose random subset of public key equations
        subset = np.random.randint(0, 2, self.n)
        
        # Compute ciphertext
        c1 = np.dot(subset, self.A) % self.q
        c2 = (np.dot(subset, self.b) + message_bit * (self.q // 2)) % self.q
        
        return c1, c2
    
    def decrypt(self, ciphertext):
        """Decrypt using secret key"""
        c1, c2 = ciphertext
        
        # Compute c2 - c1 * s
        decrypted = (c2 - np.dot(c1, self.s)) % self.q
        
        # Determine if closer to 0 or q/2
        if abs(decrypted) < abs(decrypted - self.q // 2):
            return 0
        else:
            return 1

# Demonstrate post-quantum cryptography
lwe = LWECryptosystem()
message = 1
encrypted = lwe.encrypt(message)
decrypted = lwe.decrypt(encrypted)
print(f"\nPost-quantum cryptography demo:")
print(f"Original message: {message}")
print(f"Decrypted message: {decrypted}")
print(f"Success: {message == decrypted}")

demonstrate_quantum_threat()