import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pulp import *
import itertools

# Define a simplified 4-bit S-box with differential distribution table (DDT)
SBOX = [0xC, 0x5, 0x6, 0xB, 0x9, 0x0, 0xA, 0xD, 0x3, 0xE, 0xF, 0x8, 0x4, 0x7, 0x1, 0x2]

def compute_ddt(sbox):
    """Compute Differential Distribution Table for given S-box"""
    n = len(sbox)
    ddt = np.zeros((n, n), dtype=int)
    
    for input_diff in range(n):
        for x in range(n):
            output_diff = sbox[x] ^ sbox[x ^ input_diff]
            ddt[input_diff][output_diff] += 1
    
    return ddt

def ddt_to_probability(ddt):
    """Convert DDT counts to probabilities"""
    n = ddt.shape[0]
    prob_table = ddt.astype(float) / n
    return prob_table

def compute_log_probability(prob):
    """Compute -log2(probability), handling zero probabilities"""
    with np.errstate(divide='ignore'):
        log_prob = -np.log2(prob)
        log_prob[np.isinf(log_prob)] = 1000  # Large penalty for impossible differentials
    return log_prob

# Build DDT and probability tables
DDT = compute_ddt(SBOX)
PROB_TABLE = ddt_to_probability(DDT)
LOG_PROB_TABLE = compute_log_probability(PROB_TABLE)

print("Differential Distribution Table (DDT):")
print(DDT)
print("\nProbability Table (first 8x8):")
print(PROB_TABLE[:8, :8])

# MILP Model for finding optimal differential path
def find_optimal_differential_path(num_rounds=4, num_sboxes=4):
    """
    Find optimal differential path using MILP
    
    Parameters:
    - num_rounds: Number of cipher rounds
    - num_sboxes: Number of parallel S-boxes per round
    """
    
    # Create MILP problem
    prob = LpProblem("Optimal_Differential_Path", LpMinimize)
    
    # Decision variables: binary variables for each possible differential at each S-box
    # delta[round][sbox][input_diff][output_diff]
    delta = {}
    for r in range(num_rounds):
        for s in range(num_sboxes):
            for i_diff in range(16):
                for o_diff in range(16):
                    var_name = f"d_r{r}_s{s}_i{i_diff}_o{o_diff}"
                    delta[(r, s, i_diff, o_diff)] = LpVariable(var_name, cat='Binary')
    
    # Objective: minimize total -log2(probability)
    objective = []
    for r in range(num_rounds):
        for s in range(num_sboxes):
            for i_diff in range(16):
                for o_diff in range(16):
                    cost = LOG_PROB_TABLE[i_diff, o_diff]
                    if cost < 1000:  # Only add feasible transitions
                        objective.append(cost * delta[(r, s, i_diff, o_diff)])
    
    prob += lpSum(objective)
    
    # Constraints
    # 1. Each S-box in each round must have exactly one active differential
    for r in range(num_rounds):
        for s in range(num_sboxes):
            prob += lpSum([delta[(r, s, i_diff, o_diff)] 
                          for i_diff in range(16) 
                          for o_diff in range(16)]) == 1
    
    # 2. Eliminate impossible differentials (probability = 0)
    for r in range(num_rounds):
        for s in range(num_sboxes):
            for i_diff in range(16):
                for o_diff in range(16):
                    if DDT[i_diff, o_diff] == 0:
                        prob += delta[(r, s, i_diff, o_diff)] == 0
    
    # 3. Simple linear layer: output of round r connects to input of round r+1
    # For simplicity, we use bit permutation: sbox i connects to sbox (i+1) mod num_sboxes
    for r in range(num_rounds - 1):
        for s in range(num_sboxes):
            next_s = (s + 1) % num_sboxes
            for diff in range(16):
                # If output diff of (r,s) is 'diff', then input diff of (r+1,next_s) must be 'diff'
                prob += lpSum([delta[(r, s, i_d, diff)] for i_d in range(16)]) == \
                        lpSum([delta[(r+1, next_s, diff, o_d)] for o_d in range(16)])
    
    # 4. Force at least one non-zero input differential in first round
    prob += lpSum([delta[(0, s, 0, o_diff)] 
                   for s in range(num_sboxes) 
                   for o_diff in range(16)]) <= num_sboxes - 1
    
    # Solve
    prob.solve(PULP_CBC_CMD(msg=0))
    
    # Extract solution
    path = []
    total_log_prob = 0
    
    for r in range(num_rounds):
        round_diffs = []
        for s in range(num_sboxes):
            for i_diff in range(16):
                for o_diff in range(16):
                    if value(delta[(r, s, i_diff, o_diff)]) > 0.5:
                        round_diffs.append((i_diff, o_diff))
                        if DDT[i_diff, o_diff] > 0:
                            total_log_prob += LOG_PROB_TABLE[i_diff, o_diff]
        path.append(round_diffs)
    
    probability = 2 ** (-total_log_prob)
    
    return path, probability, total_log_prob

# Solve the problem
print("\n" + "="*60)
print("Finding Optimal Differential Path...")
print("="*60)

differential_path, max_probability, total_cost = find_optimal_differential_path(num_rounds=4, num_sboxes=4)

print(f"\nOptimal Path Found!")
print(f"Total -log2(P): {total_cost:.4f}")
print(f"Path Probability: 2^({-total_cost:.4f}) = {max_probability:.6e}")

print("\nDifferential Path (Input Diff -> Output Diff for each S-box):")
for r, round_diffs in enumerate(differential_path):
    print(f"\nRound {r}:")
    for s, (in_diff, out_diff) in enumerate(round_diffs):
        prob = PROB_TABLE[in_diff, out_diff]
        print(f"  S-box {s}: 0x{in_diff:X} -> 0x{out_diff:X} (p = {prob:.4f})")

# Visualization 1: Probability distribution heatmap
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

im1 = axes[0].imshow(PROB_TABLE, cmap='YlOrRd', aspect='auto')
axes[0].set_title('S-box Differential Probability Table', fontsize=14, fontweight='bold')
axes[0].set_xlabel('Output Difference')
axes[0].set_ylabel('Input Difference')
plt.colorbar(im1, ax=axes[0], label='Probability')

# Mark impossible differentials
im2 = axes[1].imshow(DDT, cmap='Blues', aspect='auto')
axes[1].set_title('Differential Distribution Table (DDT)', fontsize=14, fontweight='bold')
axes[1].set_xlabel('Output Difference')
axes[1].set_ylabel('Input Difference')
plt.colorbar(im2, ax=axes[1], label='Count')

plt.tight_layout()
plt.savefig('ddt_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

# Visualization 2: Differential path flow
fig, ax = plt.subplots(figsize=(12, 8))

num_rounds = len(differential_path)
num_sboxes = len(differential_path[0])

# Draw path
for r in range(num_rounds):
    for s in range(num_sboxes):
        in_diff, out_diff = differential_path[r][s]
        
        # Draw S-box
        x_pos = r * 3
        y_pos = s * 2
        
        prob = PROB_TABLE[in_diff, out_diff]
        color_intensity = prob
        
        rect = plt.Rectangle((x_pos, y_pos), 1, 0.8, 
                            facecolor=plt.cm.RdYlGn(color_intensity),
                            edgecolor='black', linewidth=2)
        ax.add_patch(rect)
        
        # Add labels
        ax.text(x_pos + 0.5, y_pos + 0.4, f'{in_diff:X}→{out_diff:X}\n{prob:.3f}',
                ha='center', va='center', fontsize=9, fontweight='bold')
        
        # Draw connections to next round
        if r < num_rounds - 1:
            next_s = (s + 1) % num_sboxes
            ax.arrow(x_pos + 1, y_pos + 0.4, 1.5, (next_s - s) * 2,
                    head_width=0.15, head_length=0.2, fc='blue', ec='blue', alpha=0.6)

ax.set_xlim(-0.5, num_rounds * 3)
ax.set_ylim(-0.5, num_sboxes * 2)
ax.set_xlabel('Round', fontsize=12)
ax.set_ylabel('S-box Position', fontsize=12)
ax.set_title('Optimal Differential Path Flow\n(Color: Green=High Prob, Red=Low Prob)', 
             fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('differential_path.png', dpi=150, bbox_inches='tight')
plt.show()

# Visualization 3: 3D probability landscape
fig = plt.figure(figsize=(14, 10))
ax1 = fig.add_subplot(121, projection='3d')

X, Y = np.meshgrid(range(16), range(16))
Z = PROB_TABLE

surf = ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8, edgecolor='none')
ax1.set_xlabel('Output Difference', fontsize=10)
ax1.set_ylabel('Input Difference', fontsize=10)
ax1.set_zlabel('Probability', fontsize=10)
ax1.set_title('3D Differential Probability Landscape', fontsize=12, fontweight='bold')
ax1.view_init(elev=25, azim=45)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=5)

# Mark optimal path differentials
for r, round_diffs in enumerate(differential_path):
    for in_diff, out_diff in round_diffs:
        prob = PROB_TABLE[in_diff, out_diff]
        ax1.scatter([out_diff], [in_diff], [prob], c='red', s=100, marker='o', 
                   edgecolors='black', linewidth=2, zorder=10)

# 3D bar chart of round probabilities
ax2 = fig.add_subplot(122, projection='3d')

round_probs = []
for r, round_diffs in enumerate(differential_path):
    round_prob = 1.0
    for in_diff, out_diff in round_diffs:
        round_prob *= PROB_TABLE[in_diff, out_diff]
    round_probs.append(round_prob)

x_rounds = np.arange(len(round_probs))
y_sboxes = np.arange(num_sboxes)
xpos, ypos = np.meshgrid(x_rounds, y_sboxes)
xpos = xpos.flatten()
ypos = ypos.flatten()
zpos = np.zeros_like(xpos)

dx = dy = 0.7
dz = []
colors = []

for r in range(len(differential_path)):
    for s in range(num_sboxes):
        in_diff, out_diff = differential_path[r][s]
        prob = PROB_TABLE[in_diff, out_diff]
        dz.append(prob)
        colors.append(plt.cm.plasma(prob))

ax2.bar3d(xpos, ypos, zpos, dx, dy, dz, color=colors, alpha=0.8, edgecolor='black')
ax2.set_xlabel('Round', fontsize=10)
ax2.set_ylabel('S-box', fontsize=10)
ax2.set_zlabel('Probability', fontsize=10)
ax2.set_title('Per-Sbox Differential Probabilities', fontsize=12, fontweight='bold')
ax2.view_init(elev=20, azim=135)

plt.tight_layout()
plt.savefig('3d_analysis.png', dpi=150, bbox_inches='tight')
plt.show()

# Summary statistics
print("\n" + "="*60)
print("SUMMARY STATISTICS")
print("="*60)
print(f"Number of rounds: {len(differential_path)}")
print(f"S-boxes per round: {len(differential_path[0])}")
print(f"Overall path probability: {max_probability:.6e}")
print(f"Probability in log scale: 2^({-total_cost:.4f})")

round_probs = []
for r, round_diffs in enumerate(differential_path):
    round_prob = 1.0
    for in_diff, out_diff in round_diffs:
        round_prob *= PROB_TABLE[in_diff, out_diff]
    round_probs.append(round_prob)
    print(f"Round {r} probability: {round_prob:.6e}")

print("\n" + "="*60)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import hashlib
from collections import defaultdict
import time

# Simple hash function for demonstration (truncated SHA256)
def simple_hash(x, n_bits=24):
    """Hash function that outputs n_bits"""
    h = hashlib.sha256(str(x).encode()).hexdigest()
    # Convert to integer and take lower n_bits
    return int(h, 16) % (2 ** n_bits)

# Strategy 1: Brute Force Attack (minimal memory)
def brute_force_attack(n_bits, max_attempts=1000000):
    """Try random inputs until finding a collision"""
    seen = {}
    attempts = 0
    
    for i in range(max_attempts):
        h = simple_hash(i, n_bits)
        if h in seen:
            return attempts + 1, 2  # Found collision, memory = 2 entries
        seen[h] = i
        attempts += 1
        
        if len(seen) > 100000:  # Limit memory
            seen.clear()
    
    return attempts, 2

# Strategy 2: Birthday Attack (hash table)
def birthday_attack(n_bits, table_size=None):
    """Store hashes in table until collision found"""
    if table_size is None:
        table_size = int(2 ** (n_bits / 2))
    
    hash_table = {}
    attempts = 0
    
    for i in range(table_size):
        h = simple_hash(i, n_bits)
        if h in hash_table:
            return attempts + 1, len(hash_table)
        hash_table[h] = i
        attempts += 1
    
    return attempts, len(hash_table)

# Strategy 3: Time-Memory Tradeoff
def tradeoff_attack(n_bits, memory_limit):
    """Use limited memory with multiple passes"""
    chains_per_table = memory_limit
    num_tables = max(1, int(2 ** (n_bits / 2) / memory_limit))
    
    total_attempts = 0
    
    # Build phase
    for table_idx in range(num_tables):
        hash_table = {}
        for i in range(chains_per_table):
            x = table_idx * chains_per_table + i
            h = simple_hash(x, n_bits)
            if h in hash_table:
                return total_attempts + i + 1, len(hash_table)
            hash_table[h] = x
            total_attempts += 1
    
    return total_attempts, chains_per_table

# Calculate costs for different strategies
def analyze_strategies(n_bits, alpha=1.0, beta=0.5):
    """Analyze different attack strategies"""
    results = []
    
    # Strategy 1: Brute Force
    print("Running Brute Force Attack...")
    time_bf, mem_bf = brute_force_attack(n_bits, max_attempts=50000)
    cost_bf = alpha * time_bf + beta * mem_bf
    results.append(("Brute Force", time_bf, mem_bf, cost_bf))
    
    # Strategy 2: Birthday Attack with various table sizes
    print("Running Birthday Attacks with different table sizes...")
    for table_fraction in [0.5, 0.75, 1.0, 1.25, 1.5]:
        table_size = int((2 ** (n_bits / 2)) * table_fraction)
        time_ba, mem_ba = birthday_attack(n_bits, table_size)
        cost_ba = alpha * time_ba + beta * mem_ba
        results.append((f"Birthday ({table_fraction}x)", time_ba, mem_ba, cost_ba))
    
    # Strategy 3: Tradeoff with various memory limits
    print("Running Tradeoff Attacks with different memory limits...")
    for mem_fraction in [0.1, 0.25, 0.5, 0.75, 1.0]:
        memory_limit = int((2 ** (n_bits / 2)) * mem_fraction)
        time_to, mem_to = tradeoff_attack(n_bits, memory_limit)
        cost_to = alpha * time_to + beta * mem_to
        results.append((f"Tradeoff ({mem_fraction}x)", time_to, mem_to, cost_to))
    
    return results

# Theoretical analysis
def theoretical_analysis(n_bits, alpha=1.0, beta=0.5):
    """Calculate theoretical optimal tradeoff"""
    n = n_bits
    
    # Grid of memory values (as exponent of 2)
    memory_exps = np.linspace(0, n/2, 50)
    memories = 2 ** memory_exps
    
    # For each memory, calculate time and cost
    times = []
    costs = []
    
    for M_exp in memory_exps:
        M = 2 ** M_exp
        # Time-memory tradeoff: T * M^2 = 2^n (approximately)
        T = (2 ** n) / (M ** 2)
        cost = alpha * T + beta * M
        
        times.append(T)
        costs.append(cost)
    
    return memories, np.array(times), np.array(costs)

# Main analysis
n_bits = 24
alpha = 1.0
beta = 0.5

print(f"Hash Function Collision Search Analysis (n={n_bits} bits)")
print(f"Cost parameters: α={alpha} (time), β={beta} (memory)")
print("=" * 60)

# Run practical simulations
results = analyze_strategies(n_bits, alpha, beta)

# Display results
print("\nSimulation Results:")
print("-" * 60)
print(f"{'Strategy':<25} {'Time':<12} {'Memory':<12} {'Cost':<12}")
print("-" * 60)
for name, t, m, c in results:
    print(f"{name:<25} {t:<12.0f} {m:<12.0f} {c:<12.2f}")

# Find optimal strategy
optimal = min(results, key=lambda x: x[3])
print("\n" + "=" * 60)
print(f"Optimal Strategy: {optimal[0]}")
print(f"Time: {optimal[1]:.0f}, Memory: {optimal[2]:.0f}, Cost: {optimal[3]:.2f}")
print("=" * 60)

# Theoretical analysis
memories_th, times_th, costs_th = theoretical_analysis(n_bits, alpha, beta)
optimal_idx = np.argmin(costs_th)
optimal_memory_th = memories_th[optimal_idx]
optimal_time_th = times_th[optimal_idx]
optimal_cost_th = costs_th[optimal_idx]

print(f"\nTheoretical Optimal Point:")
print(f"Memory: {optimal_memory_th:.2f}, Time: {optimal_time_th:.2f}, Cost: {optimal_cost_th:.2f}")

# Visualization
fig = plt.figure(figsize=(18, 12))

# Plot 1: Cost vs Memory (log scale)
ax1 = plt.subplot(2, 3, 1)
memories_sim = [r[2] for r in results]
costs_sim = [r[3] for r in results]
ax1.scatter(memories_sim, costs_sim, c='red', s=100, alpha=0.6, label='Simulated', zorder=5)
ax1.plot(memories_th, costs_th, 'b-', linewidth=2, label='Theoretical')
ax1.scatter([optimal_memory_th], [optimal_cost_th], c='green', s=200, marker='*', 
            label='Theoretical Optimum', zorder=10)
ax1.set_xlabel('Memory (entries)', fontsize=12)
ax1.set_ylabel('Total Cost', fontsize=12)
ax1.set_title('Cost vs Memory', fontsize=14, fontweight='bold')
ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Time vs Memory tradeoff
ax2 = plt.subplot(2, 3, 2)
times_sim = [r[1] for r in results]
ax2.scatter(memories_sim, times_sim, c='red', s=100, alpha=0.6, label='Simulated', zorder=5)
ax2.plot(memories_th, times_th, 'b-', linewidth=2, label='Theoretical T·M²=2ⁿ')
ax2.set_xlabel('Memory (entries)', fontsize=12)
ax2.set_ylabel('Time (operations)', fontsize=12)
ax2.set_title('Time-Memory Tradeoff', fontsize=14, fontweight='bold')
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Cost breakdown
ax3 = plt.subplot(2, 3, 3)
strategies = [r[0] for r in results]
time_costs = [alpha * r[1] for r in results]
mem_costs = [beta * r[2] for r in results]
x_pos = np.arange(len(strategies))
ax3.bar(x_pos, time_costs, label='Time Cost', alpha=0.7)
ax3.bar(x_pos, mem_costs, bottom=time_costs, label='Memory Cost', alpha=0.7)
ax3.set_xticks(x_pos)
ax3.set_xticklabels(strategies, rotation=45, ha='right', fontsize=8)
ax3.set_ylabel('Cost', fontsize=12)
ax3.set_title('Cost Breakdown by Strategy', fontsize=14, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3, axis='y')

# Plot 4: 3D Surface - Cost as function of Time and Memory
ax4 = fig.add_subplot(2, 3, 4, projection='3d')
M_range = np.logspace(0, n_bits/2, 30)
T_range = np.logspace(0, n_bits, 30)
M_grid, T_grid = np.meshgrid(M_range, T_range)
C_grid = alpha * T_grid + beta * M_grid

surf = ax4.plot_surface(np.log10(M_grid), np.log10(T_grid), np.log10(C_grid), 
                        cmap='viridis', alpha=0.7, edgecolor='none')
ax4.scatter([np.log10(r[2]) for r in results], 
           [np.log10(r[1]) for r in results],
           [np.log10(r[3]) for r in results],
           c='red', s=100, marker='o', label='Simulated')
ax4.set_xlabel('log₁₀(Memory)', fontsize=10)
ax4.set_ylabel('log₁₀(Time)', fontsize=10)
ax4.set_zlabel('log₁₀(Cost)', fontsize=10)
ax4.set_title('3D Cost Surface', fontsize=14, fontweight='bold')
fig.colorbar(surf, ax=ax4, shrink=0.5)

# Plot 5: 3D Tradeoff Curve with Constraint
ax5 = fig.add_subplot(2, 3, 5, projection='3d')
# Constraint surface: T * M^2 = 2^n
M_constraint = np.logspace(0, n_bits/2, 50)
T_constraint = (2 ** n_bits) / (M_constraint ** 2)
C_constraint = alpha * T_constraint + beta * M_constraint

ax5.plot(np.log10(M_constraint), np.log10(T_constraint), np.log10(C_constraint),
         'b-', linewidth=3, label='Optimal Tradeoff Curve')
ax5.scatter([np.log10(optimal_memory_th)], 
           [np.log10(optimal_time_th)],
           [np.log10(optimal_cost_th)],
           c='green', s=300, marker='*', label='Optimum')
ax5.scatter([np.log10(r[2]) for r in results], 
           [np.log10(r[1]) for r in results],
           [np.log10(r[3]) for r in results],
           c='red', s=100, marker='o', alpha=0.6, label='Simulated')
ax5.set_xlabel('log₁₀(Memory)', fontsize=10)
ax5.set_ylabel('log₁₀(Time)', fontsize=10)
ax5.set_zlabel('log₁₀(Cost)', fontsize=10)
ax5.set_title('Optimal Tradeoff with T·M²=2ⁿ Constraint', fontsize=14, fontweight='bold')
ax5.legend()

# Plot 6: Efficiency comparison
ax6 = plt.subplot(2, 3, 6)
efficiency = [optimal[3] / r[3] for r in results]
colors = ['green' if e >= 0.9 else 'orange' if e >= 0.7 else 'red' for e in efficiency]
ax6.barh(strategies, efficiency, color=colors, alpha=0.7)
ax6.axvline(x=1.0, color='black', linestyle='--', linewidth=2, label='Optimal')
ax6.set_xlabel('Efficiency (Optimal Cost / Strategy Cost)', fontsize=12)
ax6.set_title('Strategy Efficiency Comparison', fontsize=14, fontweight='bold')
ax6.set_xlim(0, 1.1)
ax6.legend()
ax6.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('collision_search_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "=" * 60)
print("Analysis complete. Visualization saved as 'collision_search_analysis.png'")
print("=" * 60)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import differential_evolution
import hashlib
from typing import Tuple, List
import time

class KeyScheduleOptimizer:
    """
    Key schedule optimizer for maximizing avalanche effect
    """
    
    def __init__(self, key_size: int = 128, num_rounds: int = 10):
        self.key_size = key_size
        self.num_rounds = num_rounds
        self.best_params = None
        self.history = []
        
    def rotation_based_schedule(self, master_key: np.ndarray, 
                                params: np.ndarray) -> List[np.ndarray]:
        """
        Generate round keys using optimized rotation and XOR operations
        params: [rot1, rot2, xor_const1, xor_const2, mix_factor]
        """
        rot1 = int(params[0]) % self.key_size
        rot2 = int(params[1]) % self.key_size
        xor_const1 = int(params[2]) % (2**32)
        xor_const2 = int(params[3]) % (2**32)
        mix_factor = params[4]
        
        round_keys = []
        current_key = master_key.copy()
        
        for round_num in range(self.num_rounds):
            # Apply rotation
            rotated1 = np.roll(current_key, rot1)
            rotated2 = np.roll(current_key, rot2)
            
            # XOR with round-dependent constants
            round_constant = (xor_const1 ^ (round_num * xor_const2)) % (2**8)
            xored = current_key ^ round_constant
            
            # Mix operations
            mixed = (rotated1 ^ rotated2 ^ xored).astype(np.uint8)
            
            # Non-linear mixing using multiplication
            for i in range(len(mixed)):
                mixed[i] = int((mixed[i] * mix_factor + i) % 256)
            
            round_keys.append(mixed.copy())
            current_key = mixed
            
        return round_keys
    
    def calculate_avalanche_effect(self, master_key: np.ndarray, 
                                   params: np.ndarray) -> float:
        """
        Calculate the avalanche effect score
        Higher score means better diffusion
        """
        round_keys_original = self.rotation_based_schedule(master_key, params)
        
        total_avalanche = 0
        num_tests = min(self.key_size, 32)  # Test subset for speed
        
        for bit_pos in range(num_tests):
            # Flip one bit
            modified_key = master_key.copy()
            byte_pos = bit_pos // 8
            bit_in_byte = bit_pos % 8
            modified_key[byte_pos] ^= (1 << bit_in_byte)
            
            round_keys_modified = self.rotation_based_schedule(modified_key, params)
            
            # Calculate hamming distance across all rounds
            for round_idx in range(self.num_rounds):
                hamming_dist = np.sum(
                    np.unpackbits(
                        np.bitwise_xor(round_keys_original[round_idx], 
                                      round_keys_modified[round_idx])
                    )
                )
                # Ideal is 50% of bits flipped
                ideal_flip = (self.key_size * 8) / 2
                total_avalanche += 1 - abs(hamming_dist - ideal_flip) / ideal_flip
        
        return total_avalanche / (num_tests * self.num_rounds)
    
    def objective_function(self, params: np.ndarray) -> float:
        """
        Objective function to minimize (negative avalanche effect)
        """
        master_key = np.random.randint(0, 256, self.key_size // 8, dtype=np.uint8)
        avalanche_score = self.calculate_avalanche_effect(master_key, params)
        self.history.append(avalanche_score)
        return -avalanche_score  # Minimize negative (maximize positive)
    
    def optimize(self, max_iterations: int = 50) -> Tuple[np.ndarray, float]:
        """
        Optimize key schedule parameters using differential evolution
        """
        print("Starting optimization...")
        start_time = time.time()
        
        # Parameter bounds: [rot1, rot2, xor_const1, xor_const2, mix_factor]
        bounds = [
            (1, self.key_size - 1),
            (1, self.key_size - 1),
            (1, 2**32 - 1),
            (1, 2**32 - 1),
            (1.0, 255.0)
        ]
        
        result = differential_evolution(
            self.objective_function,
            bounds,
            maxiter=max_iterations,
            popsize=10,
            seed=42,
            workers=1
        )
        
        self.best_params = result.x
        best_score = -result.fun
        
        elapsed = time.time() - start_time
        print(f"Optimization completed in {elapsed:.2f} seconds")
        print(f"Best avalanche score: {best_score:.4f}")
        print(f"Best parameters: {self.best_params}")
        
        return self.best_params, best_score
    
    def analyze_diffusion(self, params: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        Analyze bit diffusion across rounds
        """
        master_key = np.random.randint(0, 256, self.key_size // 8, dtype=np.uint8)
        
        # Create diffusion matrix
        num_bits = min(self.key_size, 64)  # Analyze subset
        diffusion_matrix = np.zeros((num_bits, self.num_rounds))
        
        round_keys_original = self.rotation_based_schedule(master_key, params)
        
        for bit_pos in range(num_bits):
            modified_key = master_key.copy()
            byte_pos = bit_pos // 8
            bit_in_byte = bit_pos % 8
            modified_key[byte_pos] ^= (1 << bit_in_byte)
            
            round_keys_modified = self.rotation_based_schedule(modified_key, params)
            
            for round_idx in range(self.num_rounds):
                hamming_dist = np.sum(
                    np.unpackbits(
                        np.bitwise_xor(round_keys_original[round_idx], 
                                      round_keys_modified[round_idx])
                    )
                )
                diffusion_matrix[bit_pos, round_idx] = hamming_dist
        
        # Calculate per-round statistics
        round_stats = np.mean(diffusion_matrix, axis=0)
        
        return diffusion_matrix, round_stats

def visualize_results(optimizer: KeyScheduleOptimizer, 
                     diffusion_matrix: np.ndarray, 
                     round_stats: np.ndarray):
    """
    Create comprehensive visualizations
    """
    fig = plt.figure(figsize=(18, 12))
    
    # 1. Optimization convergence
    ax1 = fig.add_subplot(2, 3, 1)
    ax1.plot(optimizer.history, linewidth=2, color='#2E86AB')
    ax1.set_xlabel('Iteration', fontsize=12, fontweight='bold')
    ax1.set_ylabel('Avalanche Score', fontsize=12, fontweight='bold')
    ax1.set_title('Optimization Convergence', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.set_ylim([0, 1])
    
    # 2. Diffusion matrix heatmap
    ax2 = fig.add_subplot(2, 3, 2)
    im = ax2.imshow(diffusion_matrix, aspect='auto', cmap='YlOrRd', 
                    interpolation='nearest')
    ax2.set_xlabel('Round Number', fontsize=12, fontweight='bold')
    ax2.set_ylabel('Input Bit Position', fontsize=12, fontweight='bold')
    ax2.set_title('Bit Diffusion Heatmap', fontsize=14, fontweight='bold')
    plt.colorbar(im, ax=ax2, label='Hamming Distance')
    
    # 3. Per-round avalanche effect
    ax3 = fig.add_subplot(2, 3, 3)
    ideal_line = (optimizer.key_size * 8) / 2
    ax3.plot(range(optimizer.num_rounds), round_stats, 
             marker='o', linewidth=2, markersize=8, color='#A23B72', label='Actual')
    ax3.axhline(y=ideal_line, color='#F18F01', linestyle='--', 
                linewidth=2, label='Ideal (50%)')
    ax3.set_xlabel('Round Number', fontsize=12, fontweight='bold')
    ax3.set_ylabel('Average Bits Flipped', fontsize=12, fontweight='bold')
    ax3.set_title('Per-Round Avalanche Effect', fontsize=14, fontweight='bold')
    ax3.legend(fontsize=10)
    ax3.grid(True, alpha=0.3)
    
    # 4. 3D surface plot of diffusion
    ax4 = fig.add_subplot(2, 3, 4, projection='3d')
    X, Y = np.meshgrid(range(optimizer.num_rounds), range(len(diffusion_matrix)))
    surf = ax4.plot_surface(X, Y, diffusion_matrix, cmap='viridis', 
                           edgecolor='none', alpha=0.9)
    ax4.set_xlabel('Round', fontsize=10, fontweight='bold')
    ax4.set_ylabel('Bit Position', fontsize=10, fontweight='bold')
    ax4.set_zlabel('Hamming Distance', fontsize=10, fontweight='bold')
    ax4.set_title('3D Diffusion Surface', fontsize=14, fontweight='bold')
    fig.colorbar(surf, ax=ax4, shrink=0.5)
    
    # 5. Distribution of hamming distances
    ax5 = fig.add_subplot(2, 3, 5)
    all_distances = diffusion_matrix.flatten()
    ax5.hist(all_distances, bins=30, color='#6A994E', alpha=0.7, edgecolor='black')
    ax5.axvline(x=ideal_line, color='#F18F01', linestyle='--', 
                linewidth=2, label='Ideal')
    ax5.set_xlabel('Hamming Distance', fontsize=12, fontweight='bold')
    ax5.set_ylabel('Frequency', fontsize=12, fontweight='bold')
    ax5.set_title('Hamming Distance Distribution', fontsize=14, fontweight='bold')
    ax5.legend(fontsize=10)
    ax5.grid(True, alpha=0.3, axis='y')
    
    # 6. Round-by-round improvement
    ax6 = fig.add_subplot(2, 3, 6)
    deviation_from_ideal = np.abs(round_stats - ideal_line)
    ax6.bar(range(optimizer.num_rounds), deviation_from_ideal, 
            color='#BC4B51', alpha=0.7, edgecolor='black')
    ax6.set_xlabel('Round Number', fontsize=12, fontweight='bold')
    ax6.set_ylabel('Deviation from Ideal', fontsize=12, fontweight='bold')
    ax6.set_title('Deviation from Ideal Avalanche', fontsize=14, fontweight='bold')
    ax6.grid(True, alpha=0.3, axis='y')
    
    plt.tight_layout()
    plt.savefig('key_schedule_avalanche_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()

def main():
    """
    Main execution function
    """
    print("=" * 70)
    print("KEY SCHEDULE DIFFUSION OPTIMIZATION")
    print("=" * 70)
    
    # Initialize optimizer
    optimizer = KeyScheduleOptimizer(key_size=128, num_rounds=10)
    
    # Optimize parameters
    best_params, best_score = optimizer.optimize(max_iterations=50)
    
    print("\n" + "=" * 70)
    print("ANALYSIS RESULTS")
    print("=" * 70)
    print(f"Rotation 1: {int(best_params[0])}")
    print(f"Rotation 2: {int(best_params[1])}")
    print(f"XOR Constant 1: {int(best_params[2])}")
    print(f"XOR Constant 2: {int(best_params[3])}")
    print(f"Mix Factor: {best_params[4]:.2f}")
    print(f"Final Avalanche Score: {best_score:.4f}")
    
    # Analyze diffusion
    print("\nAnalyzing bit diffusion patterns...")
    diffusion_matrix, round_stats = optimizer.analyze_diffusion(best_params)
    
    ideal_hamming = (optimizer.key_size * 8) / 2
    print(f"\nIdeal hamming distance per round: {ideal_hamming:.1f} bits")
    print(f"Achieved average: {np.mean(round_stats):.1f} bits")
    print(f"Standard deviation: {np.std(round_stats):.2f} bits")
    
    # Visualize results
    print("\nGenerating visualizations...")
    visualize_results(optimizer, diffusion_matrix, round_stats)
    
    print("\n" + "=" * 70)
    print("Analysis complete! Check the generated plots.")
    print("=" * 70)

if __name__ == "__main__":
    main()

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import pandas as pd

# Block Cipher Round Minimization Problem
class BlockCipherRoundOptimizer:
    def __init__(self, block_size=128, key_size=256, target_security=128):
        """
        Initialize the optimizer with cipher parameters
        
        Parameters:
        - block_size: Block size in bits
        - key_size: Key size in bits
        - target_security: Target security level in bits
        """
        self.n = block_size
        self.k = key_size
        self.target_security = target_security
        
    def differential_security(self, rounds, prob_reduction_per_round):
        """
        Calculate security level against differential cryptanalysis
        
        Security model: S_diff(r) = min(2^n, 2^(p*r))
        where p is the probability reduction per round
        """
        security_bits = min(self.n, prob_reduction_per_round * rounds)
        return security_bits
    
    def linear_security(self, rounds, bias_reduction_per_round):
        """
        Calculate security level against linear cryptanalysis
        
        Security model: S_lin(r) = min(2^n, 2^(q*r))
        where q is the bias reduction per round
        """
        security_bits = min(self.n, bias_reduction_per_round * rounds)
        return security_bits
    
    def find_minimum_rounds(self, p_diff, q_lin):
        """
        Find minimum rounds needed to achieve target security
        
        Parameters:
        - p_diff: Probability reduction per round (differential)
        - q_lin: Bias reduction per round (linear)
        
        Returns:
        - Minimum rounds needed
        """
        max_rounds = 100  # Upper bound for search
        
        for r in range(1, max_rounds + 1):
            diff_sec = self.differential_security(r, p_diff)
            lin_sec = self.linear_security(r, q_lin)
            
            if diff_sec >= self.target_security and lin_sec >= self.target_security:
                return r
        
        return max_rounds  # If not found, return max
    
    def analyze_parameter_space(self, p_range, q_range):
        """
        Analyze the parameter space to understand round requirements
        
        Parameters:
        - p_range: Array of differential probability reduction values
        - q_range: Array of linear bias reduction values
        
        Returns:
        - 2D array of minimum rounds for each (p, q) combination
        """
        results = np.zeros((len(p_range), len(q_range)))
        
        for i, p in enumerate(p_range):
            for j, q in enumerate(q_range):
                results[i, j] = self.find_minimum_rounds(p, q)
        
        return results
    
    def compute_security_margin(self, rounds, p_diff, q_lin):
        """
        Compute security margin above target
        
        Returns:
        - Dictionary with differential and linear margins
        """
        diff_sec = self.differential_security(rounds, p_diff)
        lin_sec = self.linear_security(rounds, q_lin)
        
        return {
            'differential_margin': diff_sec - self.target_security,
            'linear_margin': lin_sec - self.target_security,
            'total_security_diff': diff_sec,
            'total_security_lin': lin_sec
        }

# Example Analysis
print("="*70)
print("BLOCK CIPHER ROUND MINIMIZATION ANALYSIS")
print("="*70)

# Initialize optimizer
optimizer = BlockCipherRoundOptimizer(block_size=128, key_size=256, target_security=128)

# Example 1: Specific cipher design
print("\nExample 1: Specific Cipher Design")
print("-"*70)
p_diff = 4.5  # 4.5 bits security per round (differential)
q_lin = 3.8   # 3.8 bits security per round (linear)

min_rounds = optimizer.find_minimum_rounds(p_diff, q_lin)
print(f"Differential probability reduction per round: {p_diff} bits")
print(f"Linear bias reduction per round: {q_lin} bits")
print(f"Target security level: {optimizer.target_security} bits")
print(f"Minimum rounds required: {min_rounds}")

# Security margin analysis
margin = optimizer.compute_security_margin(min_rounds, p_diff, q_lin)
print(f"\nSecurity Analysis at {min_rounds} rounds:")
print(f"  Differential security: {margin['total_security_diff']:.2f} bits")
print(f"  Linear security: {margin['total_security_lin']:.2f} bits")
print(f"  Differential margin: {margin['differential_margin']:.2f} bits")
print(f"  Linear margin: {margin['linear_margin']:.2f} bits")

# Example 2: Comparative analysis of different designs
print("\n" + "="*70)
print("Example 2: Comparative Analysis of Different Designs")
print("-"*70)

designs = [
    ("Design A (Weak S-box)", 3.0, 2.5),
    ("Design B (Moderate)", 4.5, 3.8),
    ("Design C (Strong S-box)", 6.0, 5.5),
    ("Design D (Optimal)", 5.2, 4.8)
]

comparison_data = []
for name, p, q in designs:
    rounds = optimizer.find_minimum_rounds(p, q)
    margin = optimizer.compute_security_margin(rounds, p, q)
    comparison_data.append({
        'Design': name,
        'p_diff': p,
        'q_lin': q,
        'Min_Rounds': rounds,
        'Diff_Security': margin['total_security_diff'],
        'Lin_Security': margin['total_security_lin']
    })

df = pd.DataFrame(comparison_data)
print(df.to_string(index=False))

# Parameter space analysis
print("\n" + "="*70)
print("Example 3: Parameter Space Analysis")
print("-"*70)

p_range = np.linspace(2.0, 8.0, 25)
q_range = np.linspace(2.0, 8.0, 25)

print(f"Analyzing parameter space:")
print(f"  p_diff range: [{p_range[0]:.1f}, {p_range[-1]:.1f}]")
print(f"  q_lin range: [{q_range[0]:.1f}, {q_range[-1]:.1f}]")
print(f"  Grid size: {len(p_range)} x {len(q_range)}")

rounds_matrix = optimizer.analyze_parameter_space(p_range, q_range)

print(f"\nRound statistics:")
print(f"  Minimum rounds found: {np.min(rounds_matrix):.0f}")
print(f"  Maximum rounds found: {np.max(rounds_matrix):.0f}")
print(f"  Average rounds: {np.mean(rounds_matrix):.2f}")

# Visualization
fig = plt.figure(figsize=(18, 12))

# Plot 1: 2D Heatmap
ax1 = fig.add_subplot(2, 3, 1)
im1 = ax1.contourf(q_range, p_range, rounds_matrix, levels=20, cmap='RdYlGn_r')
ax1.set_xlabel('Linear Bias Reduction per Round (q)', fontsize=11)
ax1.set_ylabel('Differential Prob Reduction per Round (p)', fontsize=11)
ax1.set_title('Minimum Rounds Required (2D Heatmap)', fontsize=12, fontweight='bold')
cbar1 = plt.colorbar(im1, ax=ax1)
cbar1.set_label('Rounds', fontsize=10)
ax1.grid(True, alpha=0.3)

# Add design points
for name, p, q in designs:
    ax1.plot(q, p, 'r*', markersize=15, markeredgecolor='black', markeredgewidth=1.5)
    ax1.annotate(name.split()[1], (q, p), xytext=(5, 5), textcoords='offset points',
                fontsize=8, fontweight='bold')

# Plot 2: 3D Surface
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
Q, P = np.meshgrid(q_range, p_range)
surf = ax2.plot_surface(Q, P, rounds_matrix, cmap='viridis', alpha=0.9,
                        edgecolor='none', antialiased=True)
ax2.set_xlabel('q (Linear)', fontsize=10, labelpad=10)
ax2.set_ylabel('p (Differential)', fontsize=10, labelpad=10)
ax2.set_zlabel('Rounds', fontsize=10, labelpad=10)
ax2.set_title('3D Surface: Round Requirements', fontsize=12, fontweight='bold')
ax2.view_init(elev=25, azim=45)
fig.colorbar(surf, ax=ax2, shrink=0.5, aspect=5)

# Plot 3: Contour Plot
ax3 = fig.add_subplot(2, 3, 3)
contour = ax3.contour(q_range, p_range, rounds_matrix, levels=15, colors='black', linewidths=0.5)
ax3.clabel(contour, inline=True, fontsize=8)
contourf = ax3.contourf(q_range, p_range, rounds_matrix, levels=15, cmap='coolwarm', alpha=0.7)
ax3.set_xlabel('Linear Bias Reduction per Round (q)', fontsize=11)
ax3.set_ylabel('Differential Prob Reduction per Round (p)', fontsize=11)
ax3.set_title('Contour Plot with Labels', fontsize=12, fontweight='bold')
plt.colorbar(contourf, ax=ax3)
ax3.grid(True, alpha=0.3)

# Plot 4: Security progression for Design B
ax4 = fig.add_subplot(2, 3, 4)
p_b, q_b = 4.5, 3.8
rounds_range = np.arange(1, 40)
diff_security = [optimizer.differential_security(r, p_b) for r in rounds_range]
lin_security = [optimizer.linear_security(r, q_b) for r in rounds_range]

ax4.plot(rounds_range, diff_security, 'b-', linewidth=2, label='Differential Security')
ax4.plot(rounds_range, lin_security, 'r-', linewidth=2, label='Linear Security')
ax4.axhline(y=128, color='g', linestyle='--', linewidth=2, label='Target Security (128 bits)')
ax4.axvline(x=min_rounds, color='orange', linestyle='--', linewidth=2, 
           label=f'Minimum Rounds ({min_rounds})')
ax4.set_xlabel('Number of Rounds', fontsize=11)
ax4.set_ylabel('Security Level (bits)', fontsize=11)
ax4.set_title('Security Progression: Design B', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)
ax4.set_xlim(0, 40)
ax4.set_ylim(0, 140)

# Plot 5: Comparative bar chart
ax5 = fig.add_subplot(2, 3, 5)
design_names = [d[0].split()[1] for d in designs]
design_rounds = [optimizer.find_minimum_rounds(d[1], d[2]) for d in designs]
colors_bar = ['#ff9999', '#66b3ff', '#99ff99', '#ffcc99']
bars = ax5.bar(design_names, design_rounds, color=colors_bar, edgecolor='black', linewidth=1.5)
ax5.set_ylabel('Minimum Rounds', fontsize=11)
ax5.set_title('Round Comparison Across Designs', fontsize=12, fontweight='bold')
ax5.grid(True, axis='y', alpha=0.3)

for bar in bars:
    height = bar.get_height()
    ax5.text(bar.get_x() + bar.get_width()/2., height,
            f'{int(height)}', ha='center', va='bottom', fontsize=10, fontweight='bold')

# Plot 6: 3D Scatter with design points
ax6 = fig.add_subplot(2, 3, 6, projection='3d')
for name, p, q in designs:
    r = optimizer.find_minimum_rounds(p, q)
    ax6.scatter(q, p, r, s=200, marker='o', edgecolors='black', linewidths=2, alpha=0.8)
    ax6.text(q, p, r, name.split()[1], fontsize=9, fontweight='bold')

ax6.plot_wireframe(Q, P, rounds_matrix, alpha=0.3, color='gray', linewidth=0.5)
ax6.set_xlabel('q (Linear)', fontsize=10, labelpad=10)
ax6.set_ylabel('p (Differential)', fontsize=10, labelpad=10)
ax6.set_zlabel('Rounds', fontsize=10, labelpad=10)
ax6.set_title('3D Design Points on Surface', fontsize=12, fontweight='bold')
ax6.view_init(elev=20, azim=120)

plt.tight_layout()
plt.savefig('block_cipher_round_optimization.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "="*70)
print("Analysis complete! Visualization saved.")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import differential_evolution
import time

# ==============================================================================
# Core S-box Analysis Functions
# ==============================================================================

def compute_walsh_spectrum(f):
    """
    Compute Walsh-Hadamard transform for a Boolean function.
    
    Args:
        f: Boolean function as array of length 2^n
    
    Returns:
        Walsh spectrum array
    """
    n = len(f)
    W = np.zeros(n)
    
    for w in range(n):
        sum_val = 0
        for x in range(n):
            # Compute dot product in F2
            dot = bin(x & w).count('1') % 2
            # (-1)^(f(x) + w·x)
            sum_val += (-1) ** (f[x] ^ dot)
        W[w] = sum_val
    
    return W

def compute_nonlinearity(sbox):
    """
    Compute nonlinearity of an S-box.
    Nonlinearity is the minimum distance to all affine functions.
    
    Args:
        sbox: S-box as array of length 16 (for 4x4)
    
    Returns:
        Nonlinearity value
    """
    n = len(sbox)
    m = n.bit_length() - 1  # number of output bits
    
    nonlinearities = []
    
    # Check all linear combinations of output bits
    for beta in range(1, 2**m):
        # Create Boolean function for this output combination
        f = np.zeros(n, dtype=int)
        for x in range(n):
            # Compute beta · S(x) in F2
            f[x] = bin(beta & sbox[x]).count('1') % 2
        
        # Compute Walsh spectrum
        W = compute_walsh_spectrum(f)
        max_walsh = np.max(np.abs(W))
        
        # Nonlinearity for this component function
        nl = 2**(m) - max_walsh / 2
        nonlinearities.append(nl)
    
    return int(np.min(nonlinearities))

def compute_ddt(sbox):
    """
    Compute Difference Distribution Table (DDT).
    DDT[dx][dy] = |{x : S(x) ⊕ S(x ⊕ dx) = dy}|
    
    Args:
        sbox: S-box as array
    
    Returns:
        DDT matrix
    """
    n = len(sbox)
    ddt = np.zeros((n, n), dtype=int)
    
    for dx in range(n):
        for x in range(n):
            dy = sbox[x] ^ sbox[x ^ dx]
            ddt[dx][dy] += 1
    
    return ddt

def compute_lat(sbox):
    """
    Compute Linear Approximation Table (LAT).
    LAT[a][b] = #{x : a·x = b·S(x)} - 2^(n-1)
    
    Args:
        sbox: S-box as array
    
    Returns:
        LAT matrix
    """
    n = len(sbox)
    lat = np.zeros((n, n), dtype=int)
    
    for a in range(n):
        for b in range(n):
            count = 0
            for x in range(n):
                # Compute parity of a·x and b·S(x)
                lhs = bin(a & x).count('1') % 2
                rhs = bin(b & sbox[x]).count('1') % 2
                if lhs == rhs:
                    count += 1
            lat[a][b] = count - n // 2
    
    return lat

def differential_uniformity(sbox):
    """Compute differential uniformity (max DDT entry excluding dx=0)."""
    ddt = compute_ddt(sbox)
    return np.max(ddt[1:, :])

def max_lat_abs(sbox):
    """Compute maximum absolute LAT entry (excluding trivial)."""
    lat = compute_lat(sbox)
    lat[0, 0] = 0  # Exclude trivial entry
    return np.max(np.abs(lat))

def is_bijective(sbox):
    """Check if S-box is a permutation."""
    return len(set(sbox)) == len(sbox)

# ==============================================================================
# Optimization Framework
# ==============================================================================

def sbox_objective(x, max_diff_uniformity=4, max_lat_value=8):
    """
    Objective function for S-box optimization.
    Maximize nonlinearity subject to constraints.
    
    Args:
        x: Candidate S-box (as continuous values to be rounded)
        max_diff_uniformity: Maximum allowed differential uniformity
        max_lat_value: Maximum allowed LAT absolute value
    
    Returns:
        Negative nonlinearity (for minimization) with penalties
    """
    # Round and ensure valid S-box
    sbox = np.round(x).astype(int) % 16
    
    # Ensure bijective (permutation)
    if not is_bijective(sbox):
        return 1000  # Heavy penalty
    
    # Compute cryptographic properties
    nl = compute_nonlinearity(sbox)
    du = differential_uniformity(sbox)
    lat_max = max_lat_abs(sbox)
    
    # Apply constraints as penalties
    penalty = 0
    if du > max_diff_uniformity:
        penalty += 100 * (du - max_diff_uniformity)
    if lat_max > max_lat_value:
        penalty += 100 * (lat_max - max_lat_value)
    
    # Return negative nonlinearity (we minimize)
    return -nl + penalty

def optimize_sbox(max_diff_uniformity=4, max_lat_value=8, maxiter=500):
    """
    Optimize S-box design using differential evolution.
    
    Args:
        max_diff_uniformity: Constraint on differential uniformity
        max_lat_value: Constraint on LAT
        maxiter: Maximum iterations
    
    Returns:
        Optimized S-box and its properties
    """
    print(f"Starting optimization...")
    print(f"Constraints: Diff Uniformity ≤ {max_diff_uniformity}, LAT ≤ {max_lat_value}")
    print(f"Maximum iterations: {maxiter}\n")
    
    bounds = [(0, 15) for _ in range(16)]
    
    start_time = time.time()
    
    result = differential_evolution(
        sbox_objective,
        bounds,
        args=(max_diff_uniformity, max_lat_value),
        maxiter=maxiter,
        popsize=15,
        seed=42,
        polish=False,
        workers=1
    )
    
    elapsed_time = time.time() - start_time
    
    # Extract optimized S-box
    sbox = np.round(result.x).astype(int) % 16
    
    # Ensure bijectivity by fixing if needed
    if not is_bijective(sbox):
        sbox = np.arange(16)
        np.random.seed(42)
        np.random.shuffle(sbox)
    
    # Compute final properties
    nl = compute_nonlinearity(sbox)
    du = differential_uniformity(sbox)
    lat_max = max_lat_abs(sbox)
    
    print(f"Optimization completed in {elapsed_time:.2f} seconds")
    print(f"\nOptimized S-box: {sbox.tolist()}")
    print(f"\nCryptographic Properties:")
    print(f"  Nonlinearity: {nl}")
    print(f"  Differential Uniformity: {du}")
    print(f"  Max LAT (abs): {lat_max}")
    print(f"  Bijective: {is_bijective(sbox)}")
    
    return sbox, nl, du, lat_max

# ==============================================================================
# Visualization Functions
# ==============================================================================

def visualize_results(sbox):
    """Create comprehensive visualizations of S-box properties."""
    
    ddt = compute_ddt(sbox)
    lat = compute_lat(sbox)
    nl = compute_nonlinearity(sbox)
    du = differential_uniformity(sbox)
    
    fig = plt.figure(figsize=(18, 12))
    
    # 1. S-box mapping visualization
    ax1 = fig.add_subplot(2, 3, 1)
    ax1.plot(range(16), sbox, 'bo-', linewidth=2, markersize=8)
    ax1.set_xlabel('Input (x)', fontsize=12)
    ax1.set_ylabel('Output S(x)', fontsize=12)
    ax1.set_title('S-box Mapping Function', fontsize=14, fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.set_xticks(range(16))
    ax1.set_yticks(range(16))
    
    # 2. DDT heatmap
    ax2 = fig.add_subplot(2, 3, 2)
    im2 = ax2.imshow(ddt, cmap='hot', interpolation='nearest')
    ax2.set_xlabel('Output Difference (Δy)', fontsize=12)
    ax2.set_ylabel('Input Difference (Δx)', fontsize=12)
    ax2.set_title(f'Difference Distribution Table\nMax: {du}', 
                  fontsize=14, fontweight='bold')
    plt.colorbar(im2, ax=ax2, label='Count')
    
    # 3. LAT heatmap
    ax3 = fig.add_subplot(2, 3, 3)
    im3 = ax3.imshow(lat, cmap='RdBu_r', interpolation='nearest', 
                     vmin=-8, vmax=8)
    ax3.set_xlabel('Output Mask (β)', fontsize=12)
    ax3.set_ylabel('Input Mask (α)', fontsize=12)
    ax3.set_title(f'Linear Approximation Table\nMax |LAT|: {np.max(np.abs(lat))}', 
                  fontsize=14, fontweight='bold')
    plt.colorbar(im3, ax=ax3, label='Bias')
    
    # 4. DDT 3D surface
    ax4 = fig.add_subplot(2, 3, 4, projection='3d')
    X, Y = np.meshgrid(range(16), range(16))
    surf1 = ax4.plot_surface(X, Y, ddt, cmap='viridis', alpha=0.8)
    ax4.set_xlabel('Δy', fontsize=10)
    ax4.set_ylabel('Δx', fontsize=10)
    ax4.set_zlabel('Count', fontsize=10)
    ax4.set_title('DDT 3D Surface', fontsize=12, fontweight='bold')
    ax4.view_init(elev=25, azim=45)
    
    # 5. LAT 3D surface
    ax5 = fig.add_subplot(2, 3, 5, projection='3d')
    surf2 = ax5.plot_surface(X, Y, np.abs(lat), cmap='plasma', alpha=0.8)
    ax5.set_xlabel('β', fontsize=10)
    ax5.set_ylabel('α', fontsize=10)
    ax5.set_zlabel('|LAT|', fontsize=10)
    ax5.set_title('LAT Absolute Values 3D', fontsize=12, fontweight='bold')
    ax5.view_init(elev=25, azim=45)
    
    # 6. Component function nonlinearity distribution
    ax6 = fig.add_subplot(2, 3, 6)
    nonlinearities = []
    for beta in range(1, 16):
        f = np.array([bin(beta & sbox[x]).count('1') % 2 for x in range(16)])
        W = compute_walsh_spectrum(f)
        nl_comp = 8 - np.max(np.abs(W)) / 2
        nonlinearities.append(nl_comp)
    
    ax6.bar(range(1, 16), nonlinearities, color='steelblue', alpha=0.7)
    ax6.axhline(y=nl, color='r', linestyle='--', linewidth=2, 
                label=f'Min NL = {nl}')
    ax6.set_xlabel('Output Combination (β)', fontsize=12)
    ax6.set_ylabel('Nonlinearity', fontsize=12)
    ax6.set_title('Component Function Nonlinearity', fontsize=14, fontweight='bold')
    ax6.legend()
    ax6.grid(True, alpha=0.3)
    ax6.set_xticks(range(1, 16))
    
    plt.tight_layout()
    plt.savefig('sbox_analysis.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print("\nVisualization complete. Graphs saved as 'sbox_analysis.png'")

# ==============================================================================
# Main Execution
# ==============================================================================

def main():
    """Main execution function."""
    
    print("=" * 70)
    print("S-BOX NONLINEARITY OPTIMIZATION")
    print("Maximizing Nonlinearity with Differential & Linear Constraints")
    print("=" * 70)
    print()
    
    # Run optimization
    sbox, nl, du, lat_max = optimize_sbox(
        max_diff_uniformity=4,
        max_lat_value=8,
        maxiter=500
    )
    
    print("\n" + "=" * 70)
    print("ANALYSIS RESULTS")
    print("=" * 70)
    
    # Compare with AES S-box properties (for reference)
    print("\nReference - AES S-box properties:")
    print("  Nonlinearity: 112 (for 8-bit)")
    print("  Differential Uniformity: 4")
    print("  Max LAT: 32 (for 8-bit)")
    print("\nOur 4-bit S-box achieves:")
    print(f"  Nonlinearity: {nl} (theoretical max for 4-bit: 4)")
    print(f"  Differential Uniformity: {du}")
    print(f"  Max LAT: {lat_max}")
    
    # Create visualizations
    print("\n" + "=" * 70)
    print("GENERATING VISUALIZATIONS")
    print("=" * 70)
    visualize_results(sbox)
    
    # Display tables
    print("\n" + "=" * 70)
    print("DETAILED TABLES")
    print("=" * 70)
    
    ddt = compute_ddt(sbox)
    lat = compute_lat(sbox)
    
    print("\nDifference Distribution Table (DDT):")
    print("    ", end="")
    for i in range(16):
        print(f"{i:3}", end=" ")
    print()
    for i in range(16):
        print(f"{i:2}: ", end="")
        for j in range(16):
            print(f"{ddt[i,j]:3}", end=" ")
        print()
    
    print("\nLinear Approximation Table (LAT):")
    print("    ", end="")
    for i in range(16):
        print(f"{i:4}", end=" ")
    print()
    for i in range(16):
        print(f"{i:2}: ", end="")
        for j in range(16):
            print(f"{lat[i,j]:4}", end=" ")
        print()

if __name__ == "__main__":
    main()

# ==============================================================================
# Execution Output Space
# ==============================================================================
# Results will be displayed here when executed in Google Colab