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

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

# Security estimation function
def estimate_security(n, log_q, sigma=3.2):
    """
    Estimate bit-security level using simplified LWE hardness formula
    
    Parameters:
    n: polynomial degree
    log_q: log2 of coefficient modulus
    sigma: noise standard deviation
    
    Returns:
    Estimated security level in bits
    """
    if n <= 0 or log_q <= 0:
        return 0
    return (n * log_q) / (sigma**2 * np.log2(n))

# Communication cost
def communication_cost(n, log_q):
    """Communication cost proportional to n * log_q"""
    return n * log_q

# Computation cost
def computation_cost(n):
    """Computation cost proportional to n * log(n)"""
    return n * np.log2(n)

# Optimization function
def optimize_parameters(target_security=128, sigma=3.2):
    """
    Find optimal (n, q) parameters that minimize costs under security constraint
    
    Parameters:
    target_security: target security level in bits (default: 128)
    sigma: noise standard deviation
    
    Returns:
    DataFrame with valid parameter combinations and their costs
    """
    # Candidate polynomial degrees (powers of 2)
    n_candidates = [2**i for i in range(10, 16)]  # 1024 to 32768
    
    results = []
    
    for n in n_candidates:
        # Search for minimum log_q that satisfies security constraint
        # Start from a reasonable minimum and search upward
        log_q_min = 20
        log_q_max = 1000
        
        for log_q in range(log_q_min, log_q_max):
            security = estimate_security(n, log_q, sigma)
            
            if security >= target_security:
                # Found minimum log_q for this n
                comm_cost = communication_cost(n, log_q)
                comp_cost = computation_cost(n)
                
                results.append({
                    'n': n,
                    'log_q': log_q,
                    'security': security,
                    'communication_cost': comm_cost,
                    'computation_cost': comp_cost,
                    'total_cost': comm_cost + comp_cost  # Combined cost
                })
                break
    
    df = pd.DataFrame(results)
    return df

# Generate detailed parameter space for visualization
def generate_parameter_space(n_range, log_q_range, sigma=3.2):
    """Generate full parameter space for 3D visualization"""
    N, LOG_Q = np.meshgrid(n_range, log_q_range)
    SECURITY = np.zeros_like(N, dtype=float)
    COMM_COST = np.zeros_like(N, dtype=float)
    
    for i in range(N.shape[0]):
        for j in range(N.shape[1]):
            n = N[i, j]
            log_q = LOG_Q[i, j]
            SECURITY[i, j] = estimate_security(n, log_q, sigma)
            COMM_COST[i, j] = communication_cost(n, log_q)
    
    return N, LOG_Q, SECURITY, COMM_COST

# Main execution
print("=" * 60)
print("Polynomial Modulus and Degree Parameter Optimization")
print("=" * 60)

# Set parameters
TARGET_SECURITY = 128
SIGMA = 3.2

print(f"\nTarget Security Level: {TARGET_SECURITY} bits")
print(f"Noise Standard Deviation (σ): {SIGMA}")
print("\nSearching for optimal parameters...\n")

# Optimize parameters
results_df = optimize_parameters(target_security=TARGET_SECURITY, sigma=SIGMA)

# Display results
print("Valid Parameter Combinations:")
print("=" * 80)
print(results_df.to_string(index=False))
print("=" * 80)

# Find optimal solution
optimal_idx = results_df['total_cost'].idxmin()
optimal_params = results_df.iloc[optimal_idx]

print(f"\n{'*' * 60}")
print("OPTIMAL PARAMETERS:")
print(f"{'*' * 60}")
print(f"  Polynomial Degree (n): {optimal_params['n']:.0f}")
print(f"  Log2(Modulus) (log q): {optimal_params['log_q']:.0f}")
print(f"  Achieved Security: {optimal_params['security']:.2f} bits")
print(f"  Communication Cost: {optimal_params['communication_cost']:.0f}")
print(f"  Computation Cost: {optimal_params['computation_cost']:.2f}")
print(f"  Total Cost: {optimal_params['total_cost']:.2f}")
print(f"{'*' * 60}\n")

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

# Plot 1: Parameter Trade-offs
ax1 = fig.add_subplot(2, 3, 1)
ax1.scatter(results_df['n'], results_df['log_q'], 
           s=200, c=results_df['security'], cmap='viridis', 
           edgecolors='black', linewidth=1.5)
ax1.scatter(optimal_params['n'], optimal_params['log_q'], 
           s=400, marker='*', c='red', edgecolors='black', linewidth=2,
           label='Optimal', zorder=5)
ax1.set_xlabel('Polynomial Degree (n)', fontsize=12, fontweight='bold')
ax1.set_ylabel('Log2(Modulus)', fontsize=12, fontweight='bold')
ax1.set_title('Parameter Space (Security Level)', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=10)
cbar1 = plt.colorbar(ax1.collections[0], ax=ax1)
cbar1.set_label('Security (bits)', fontsize=10)

# Plot 2: Communication Cost vs n
ax2 = fig.add_subplot(2, 3, 2)
ax2.plot(results_df['n'], results_df['communication_cost'], 
        'o-', linewidth=2, markersize=8, color='blue', label='Communication Cost')
ax2.plot(optimal_params['n'], optimal_params['communication_cost'], 
        '*', markersize=20, color='red', label='Optimal', zorder=5)
ax2.set_xlabel('Polynomial Degree (n)', fontsize=12, fontweight='bold')
ax2.set_ylabel('Communication Cost', fontsize=12, fontweight='bold')
ax2.set_title('Communication Cost vs Polynomial Degree', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend(fontsize=10)

# Plot 3: Computation Cost vs n
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(results_df['n'], results_df['computation_cost'], 
        's-', linewidth=2, markersize=8, color='green', label='Computation Cost')
ax3.plot(optimal_params['n'], optimal_params['computation_cost'], 
        '*', markersize=20, color='red', label='Optimal', zorder=5)
ax3.set_xlabel('Polynomial Degree (n)', fontsize=12, fontweight='bold')
ax3.set_ylabel('Computation Cost', fontsize=12, fontweight='bold')
ax3.set_title('Computation Cost vs Polynomial Degree', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.legend(fontsize=10)

# Plot 4: Total Cost Comparison
ax4 = fig.add_subplot(2, 3, 4)
x_pos = np.arange(len(results_df))
bars = ax4.bar(x_pos, results_df['total_cost'], 
              color=['red' if i == optimal_idx else 'skyblue' for i in range(len(results_df))],
              edgecolor='black', linewidth=1.5)
ax4.set_xlabel('Configuration', fontsize=12, fontweight='bold')
ax4.set_ylabel('Total Cost', fontsize=12, fontweight='bold')
ax4.set_title('Total Cost for Each Configuration', fontsize=14, fontweight='bold')
ax4.set_xticks(x_pos)
ax4.set_xticklabels([f"n={int(n)}" for n in results_df['n']], rotation=45)
ax4.grid(True, alpha=0.3, axis='y')

# Plot 5: Security Level Comparison
ax5 = fig.add_subplot(2, 3, 5)
ax5.barh(results_df['n'].astype(str), results_df['security'], 
        color=['red' if i == optimal_idx else 'lightcoral' for i in range(len(results_df))],
        edgecolor='black', linewidth=1.5)
ax5.axvline(x=TARGET_SECURITY, color='blue', linestyle='--', linewidth=2, label=f'Target: {TARGET_SECURITY} bits')
ax5.set_xlabel('Security Level (bits)', fontsize=12, fontweight='bold')
ax5.set_ylabel('Polynomial Degree (n)', fontsize=12, fontweight='bold')
ax5.set_title('Achieved Security Levels', fontsize=14, fontweight='bold')
ax5.grid(True, alpha=0.3, axis='x')
ax5.legend(fontsize=10)

# Plot 6: 3D Surface Plot - Security Landscape
ax6 = fig.add_subplot(2, 3, 6, projection='3d')

# Generate parameter space for visualization
n_range = np.array([2**i for i in range(10, 16)])
log_q_range = np.arange(20, 300, 5)
N, LOG_Q, SECURITY, _ = generate_parameter_space(n_range, log_q_range, SIGMA)

# Create surface plot
surf = ax6.plot_surface(np.log2(N), LOG_Q, SECURITY, 
                        cmap=cm.coolwarm, alpha=0.8, 
                        edgecolor='none', antialiased=True)

# Mark the security constraint plane
xx, yy = np.meshgrid(np.log2(n_range), log_q_range)
zz = np.ones_like(xx) * TARGET_SECURITY
ax6.plot_surface(xx, yy, zz, alpha=0.3, color='green')

# Mark optimal point
ax6.scatter([np.log2(optimal_params['n'])], 
           [optimal_params['log_q']], 
           [optimal_params['security']], 
           color='red', s=200, marker='*', 
           edgecolors='black', linewidth=2, zorder=10)

ax6.set_xlabel('Log2(n)', fontsize=10, fontweight='bold')
ax6.set_ylabel('Log2(q)', fontsize=10, fontweight='bold')
ax6.set_zlabel('Security (bits)', fontsize=10, fontweight='bold')
ax6.set_title('3D Security Landscape', fontsize=14, fontweight='bold')
fig.colorbar(surf, ax=ax6, shrink=0.5, aspect=5)

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

print("\nVisualization complete! Graph saved as 'parameter_optimization.png'")

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize_scalar
from matplotlib import cm

# Set random seed for reproducibility
np.random.seed(42)

# Problem parameters
n = 100  # Lattice dimension
q = 1009  # Modulus for LWE
sigma = 3.0  # Noise parameter

# BKZ parameters range
beta_min = 10
beta_max = 80
beta_range = np.arange(beta_min, beta_max + 1)

# Time budget scenarios (in seconds)
time_budgets = [1, 10, 100, 1000, 10000]

print("="*70)
print("BKZ Block Size Optimization Analysis")
print("="*70)
print(f"Lattice dimension (n): {n}")
print(f"Modulus (q): {q}")
print(f"Noise parameter (σ): {sigma}")
print(f"Block size range (β): [{beta_min}, {beta_max}]")
print("="*70)

def compute_root_hermite_factor(beta, n):
    """
    Compute the root Hermite factor achieved by BKZ-β
    Using the approximation: δ ≈ ((πβ)^(1/β) * β / (2πe))^(1/(2(β-1)))
    """
    if beta <= 1:
        return 1.0
    # Chen-Nguyen approximation
    delta = ((np.pi * beta) ** (1.0 / beta) * beta / (2 * np.pi * np.e)) ** (1.0 / (2 * (beta - 1)))
    return delta

def compute_computation_time(beta):
    """
    Estimate BKZ-β running time
    T(β) ≈ c * 2^(k*β) where k ≈ 0.18 for practical implementations
    """
    c = 0.001  # Base constant (seconds)
    k = 0.18   # Exponential factor
    return c * (2 ** (k * beta))

def compute_attack_success_probability(delta, n, q, sigma):
    """
    Estimate attack success probability based on root Hermite factor
    The attack succeeds if the shortest vector found is short enough
    Success probability increases as δ decreases (better reduction)
    """
    # Expected shortest vector length after BKZ reduction
    expected_length = delta ** n * (q ** (1/n))
    
    # Target length for successful attack (related to noise)
    target_length = sigma * np.sqrt(n)
    
    # Success probability using a sigmoid function
    # Higher probability when expected_length < target_length
    ratio = expected_length / target_length
    prob = 1.0 / (1.0 + np.exp(5 * (ratio - 1)))
    
    return prob

def compute_expected_success_rate(beta, time_budget, n, q, sigma):
    """
    Compute expected success rate considering time budget
    Expected success = (time_budget / T(β)) * P_success(β)
    This represents how many successful attacks we expect within the time budget
    """
    comp_time = compute_computation_time(beta)
    if comp_time > time_budget:
        # Can't even complete one attempt
        return 0.0
    
    delta = compute_root_hermite_factor(beta, n)
    prob_success = compute_attack_success_probability(delta, n, q, sigma)
    
    # Number of attempts possible within time budget
    num_attempts = time_budget / comp_time
    
    # Expected number of successes
    expected_successes = num_attempts * prob_success
    
    return expected_successes

# Compute metrics for all block sizes
deltas = np.array([compute_root_hermite_factor(beta, n) for beta in beta_range])
comp_times = np.array([compute_computation_time(beta) for beta in beta_range])
success_probs = np.array([compute_attack_success_probability(delta, n, q, sigma) 
                          for delta in deltas])

print("\nComputed Metrics Sample (every 10 block sizes):")
print(f"{'β':<5} {'δ':<12} {'Time (s)':<15} {'P(success)':<12}")
print("-" * 50)
for i in range(0, len(beta_range), 10):
    print(f"{beta_range[i]:<5} {deltas[i]:<12.6f} {comp_times[i]:<15.4e} {success_probs[i]:<12.6f}")

# Find optimal block size for each time budget
optimal_results = []
print("\n" + "="*70)
print("Optimization Results for Different Time Budgets")
print("="*70)

for time_budget in time_budgets:
    expected_rates = np.array([compute_expected_success_rate(beta, time_budget, n, q, sigma) 
                               for beta in beta_range])
    
    optimal_idx = np.argmax(expected_rates)
    optimal_beta = beta_range[optimal_idx]
    optimal_rate = expected_rates[optimal_idx]
    
    optimal_results.append({
        'time_budget': time_budget,
        'optimal_beta': optimal_beta,
        'expected_rate': optimal_rate,
        'computation_time': compute_computation_time(optimal_beta),
        'success_prob': success_probs[optimal_idx]
    })
    
    print(f"\nTime Budget: {time_budget} seconds")
    print(f"  Optimal β: {optimal_beta}")
    print(f"  Expected successes: {optimal_rate:.4f}")
    print(f"  Single run time: {compute_computation_time(optimal_beta):.4f} seconds")
    print(f"  Single run success probability: {success_probs[optimal_idx]:.6f}")
    print(f"  Number of possible runs: {time_budget / compute_computation_time(optimal_beta):.2f}")

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

# Plot 1: Computation Time vs Block Size
ax1 = fig.add_subplot(2, 3, 1)
ax1.semilogy(beta_range, comp_times, 'b-', linewidth=2, label='Computation Time')
ax1.set_xlabel('Block Size (β)', fontsize=12)
ax1.set_ylabel('Time (seconds, log scale)', fontsize=12)
ax1.set_title('BKZ Computation Time vs Block Size', fontsize=13, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Plot 2: Root Hermite Factor vs Block Size
ax2 = fig.add_subplot(2, 3, 2)
ax2.plot(beta_range, deltas, 'r-', linewidth=2, label='Root Hermite Factor δ')
ax2.set_xlabel('Block Size (β)', fontsize=12)
ax2.set_ylabel('Root Hermite Factor (δ)', fontsize=12)
ax2.set_title('Root Hermite Factor vs Block Size', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend()

# Plot 3: Success Probability vs Block Size
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(beta_range, success_probs, 'g-', linewidth=2, label='Success Probability')
ax3.set_xlabel('Block Size (β)', fontsize=12)
ax3.set_ylabel('Attack Success Probability', fontsize=12)
ax3.set_title('Attack Success Probability vs Block Size', fontsize=13, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.legend()

# Plot 4: Expected Success Rate for Different Time Budgets
ax4 = fig.add_subplot(2, 3, 4)
for time_budget in time_budgets:
    expected_rates = np.array([compute_expected_success_rate(beta, time_budget, n, q, sigma) 
                               for beta in beta_range])
    ax4.plot(beta_range, expected_rates, linewidth=2, label=f'T={time_budget}s', alpha=0.8)
ax4.set_xlabel('Block Size (β)', fontsize=12)
ax4.set_ylabel('Expected Number of Successes', fontsize=12)
ax4.set_title('Expected Success Rate vs Block Size\n(Different Time Budgets)', 
              fontsize=13, fontweight='bold')
ax4.grid(True, alpha=0.3)
ax4.legend()

# Plot 5: Optimal Block Size vs Time Budget
ax5 = fig.add_subplot(2, 3, 5)
budgets = [r['time_budget'] for r in optimal_results]
optimal_betas = [r['optimal_beta'] for r in optimal_results]
ax5.semilogx(budgets, optimal_betas, 'o-', linewidth=2, markersize=8, color='purple')
ax5.set_xlabel('Time Budget (seconds, log scale)', fontsize=12)
ax5.set_ylabel('Optimal Block Size (β)', fontsize=12)
ax5.set_title('Optimal Block Size vs Time Budget', fontsize=13, fontweight='bold')
ax5.grid(True, alpha=0.3)

# Plot 6: Trade-off Summary
ax6 = fig.add_subplot(2, 3, 6)
ax6_twin = ax6.twinx()
ax6.plot(beta_range, comp_times, 'b-', linewidth=2, label='Computation Time')
ax6_twin.plot(beta_range, success_probs, 'r-', linewidth=2, label='Success Probability')
ax6.set_xlabel('Block Size (β)', fontsize=12)
ax6.set_ylabel('Time (seconds)', fontsize=12, color='b')
ax6_twin.set_ylabel('Success Probability', fontsize=12, color='r')
ax6.set_title('Time vs Success Trade-off', fontsize=13, fontweight='bold')
ax6.tick_params(axis='y', labelcolor='b')
ax6_twin.tick_params(axis='y', labelcolor='r')
ax6.set_yscale('log')
ax6.grid(True, alpha=0.3)
lines1, labels1 = ax6.get_legend_handles_labels()
lines2, labels2 = ax6_twin.get_legend_handles_labels()
ax6.legend(lines1 + lines2, labels1 + labels2, loc='center left')

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

# 3D Visualization: Expected Success Rate Surface
fig_3d = plt.figure(figsize=(16, 6))

# 3D Plot 1: Success Rate Surface (Time Budget vs Block Size)
ax_3d1 = fig_3d.add_subplot(1, 2, 1, projection='3d')

# Create mesh for 3D surface
time_budget_range = np.logspace(0, 4, 50)  # 1 to 10000 seconds
beta_mesh, time_mesh = np.meshgrid(beta_range, time_budget_range)
success_rate_mesh = np.zeros_like(beta_mesh, dtype=float)

for i in range(len(time_budget_range)):
    for j in range(len(beta_range)):
        success_rate_mesh[i, j] = compute_expected_success_rate(
            beta_range[j], time_budget_range[i], n, q, sigma
        )

surf1 = ax_3d1.plot_surface(beta_mesh, np.log10(time_mesh), success_rate_mesh,
                            cmap=cm.viridis, alpha=0.9, edgecolor='none')
ax_3d1.set_xlabel('Block Size (β)', fontsize=11, labelpad=10)
ax_3d1.set_ylabel('log₁₀(Time Budget)', fontsize=11, labelpad=10)
ax_3d1.set_zlabel('Expected Successes', fontsize=11, labelpad=10)
ax_3d1.set_title('Expected Success Rate Surface\n(Time Budget vs Block Size)', 
                 fontsize=13, fontweight='bold', pad=20)
fig_3d.colorbar(surf1, ax=ax_3d1, shrink=0.5, aspect=5)

# Mark optimal points
for result in optimal_results:
    ax_3d1.scatter([result['optimal_beta']], 
                   [np.log10(result['time_budget'])], 
                   [result['expected_rate']],
                   color='red', s=100, marker='o', edgecolors='black', linewidths=2)

# 3D Plot 2: Trade-off Landscape (Root Hermite Factor vs Time vs Success)
ax_3d2 = fig_3d.add_subplot(1, 2, 2, projection='3d')

# Create data for scatter plot
scatter_betas = beta_range[::2]  # Sample every other point for clarity
scatter_deltas = deltas[::2]
scatter_times = comp_times[::2]
scatter_probs = success_probs[::2]

scatter = ax_3d2.scatter(scatter_deltas, np.log10(scatter_times), scatter_probs,
                        c=scatter_betas, cmap=cm.plasma, s=100, alpha=0.8,
                        edgecolors='black', linewidths=0.5)
ax_3d2.set_xlabel('Root Hermite Factor (δ)', fontsize=11, labelpad=10)
ax_3d2.set_ylabel('log₁₀(Time)', fontsize=11, labelpad=10)
ax_3d2.set_zlabel('Success Probability', fontsize=11, labelpad=10)
ax_3d2.set_title('Trade-off Landscape\n(Quality vs Time vs Success)', 
                 fontsize=13, fontweight='bold', pad=20)
cbar = fig_3d.colorbar(scatter, ax=ax_3d2, shrink=0.5, aspect=5)
cbar.set_label('Block Size (β)', fontsize=10)

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

# Summary statistics
print("\n" + "="*70)
print("Summary Statistics")
print("="*70)
print(f"Minimum computation time: {comp_times.min():.6f} seconds (β={beta_range[np.argmin(comp_times)]})")
print(f"Maximum computation time: {comp_times.max():.2f} seconds (β={beta_range[np.argmax(comp_times)]})")
print(f"Minimum root Hermite factor: {deltas.min():.6f} (β={beta_range[np.argmin(deltas)]})")
print(f"Maximum root Hermite factor: {deltas.max():.6f} (β={beta_range[np.argmax(deltas)]})")
print(f"Minimum success probability: {success_probs.min():.8f} (β={beta_range[np.argmin(success_probs)]})")
print(f"Maximum success probability: {success_probs.max():.6f} (β={beta_range[np.argmax(success_probs)]})")

print("\n" + "="*70)
print("Key Insights")
print("="*70)
print("1. Computation time grows exponentially with block size: T(β) ~ 2^(0.18β)")
print("2. Root Hermite factor decreases (improves) with larger block sizes")
print("3. Attack success probability increases with better reduction quality")
print("4. Optimal block size increases logarithmically with time budget")
print("5. The trade-off is non-trivial: maximum success probability doesn't always")
print("   yield maximum expected success rate within a time budget")
print("="*70)

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import erfc
from scipy.optimize import brentq
import warnings
warnings.filterwarnings('ignore')

# Set random seed for reproducibility
np.random.seed(42)

# LWE Parameters
n = 50  # dimension
q = 1021  # modulus (prime)
sigma = 3.2  # error standard deviation

print("=== LWE Sample Count Optimization ===")
print(f"Parameters: n={n}, q={q}, σ={sigma}")
print()

def estimate_security_bits(n, q, sigma, m):
    """
    Estimate security level in bits based on LWE parameters.
    This uses a simplified model of the primal attack.
    
    The security estimate is based on:
    - Lattice dimension: d = n + m
    - Root Hermite factor: δ = (β/d)^(1/(2β-d)) where β is BKZ block size
    - Expected shortest vector length vs required length
    """
    d = n + m  # lattice dimension
    
    # Estimate BKZ block size needed (simplified)
    # For actual attacks, this would use more sophisticated models
    log_q = np.log2(q)
    log_sigma = np.log2(sigma)
    
    # Gaussian heuristic for shortest vector
    delta = 1.0045  # root Hermite factor for moderate security
    
    # Security bits approximation
    # Based on the fact that attack cost is roughly 2^(0.292*β) for BKZ-β
    beta = min(d, int(d * 0.7))  # effective block size
    security = 0.292 * beta - 16.4 + log_sigma
    
    return max(0, security)

def attack_success_probability(n, q, sigma, m):
    """
    Model the probability of successful attack given m samples.
    
    This models the success probability based on:
    1. Information-theoretic requirements (need m > n)
    2. Statistical advantage from having more samples
    3. Error accumulation effects
    """
    if m <= n:
        # Insufficient samples - very low success probability
        return 0.01 * (m / n)
    
    # Calculate effective advantage
    # More samples = more information = higher success
    excess_samples = m - n
    
    # Model: success probability increases with sample count
    # but with diminishing returns
    log_advantage = excess_samples / (2 * sigma**2)
    
    # Apply sigmoid-like function for realistic probability
    # Success probability increases as we get more samples beyond minimum
    z = (log_advantage - 2.0) / 0.5
    prob = 1.0 / (1.0 + np.exp(-z))
    
    # Factor in dimension effects
    dim_factor = np.exp(-n / 200.0)
    prob = prob * (1.0 - dim_factor) + dim_factor * 0.5
    
    return min(0.999, max(0.001, prob))

def find_optimal_samples(n, q, sigma, target_prob):
    """
    Find the minimum number of samples needed to achieve target success probability.
    Uses binary search for efficiency.
    """
    # Start with theoretical minimum
    m_min = n + 1
    m_max = n + 500
    
    # Check if target is achievable
    if attack_success_probability(n, q, sigma, m_max) < target_prob:
        print(f"Warning: Target probability {target_prob} may not be achievable with m <= {m_max}")
        return m_max
    
    # Binary search
    while m_max - m_min > 1:
        m_mid = (m_min + m_max) // 2
        prob = attack_success_probability(n, q, sigma, m_mid)
        
        if prob < target_prob:
            m_min = m_mid
        else:
            m_max = m_mid
    
    return m_max

# Find optimal sample count for different target probabilities
target_probs = [0.50, 0.75, 0.90, 0.95, 0.99]
optimal_samples = []

print("Target Probability -> Optimal Sample Count:")
print("-" * 50)

for target_prob in target_probs:
    m_opt = find_optimal_samples(n, q, sigma, target_prob)
    optimal_samples.append(m_opt)
    actual_prob = attack_success_probability(n, q, sigma, m_opt)
    security = estimate_security_bits(n, q, sigma, m_opt)
    print(f"P_target = {target_prob:.2f} -> m = {m_opt:3d} (actual P = {actual_prob:.4f}, security ≈ {security:.1f} bits)")

print()

# Generate data for visualization
m_range = np.arange(n, n + 200, 1)
success_probs = [attack_success_probability(n, q, sigma, m) for m in m_range]
security_bits = [estimate_security_bits(n, q, sigma, m) for m in m_range]

# Create 2D plots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Success Probability vs Sample Count
ax1.plot(m_range, success_probs, 'b-', linewidth=2, label='Success Probability')
ax1.axhline(y=0.95, color='r', linestyle='--', linewidth=1.5, label='Target (0.95)')
ax1.axvline(x=optimal_samples[3], color='g', linestyle='--', linewidth=1.5, 
            label=f'Optimal m={optimal_samples[3]}')
ax1.fill_between(m_range, 0, success_probs, alpha=0.3)
ax1.set_xlabel('Number of Samples (m)', fontsize=12)
ax1.set_ylabel('Attack Success Probability', fontsize=12)
ax1.set_title('LWE Attack Success Probability vs Sample Count', fontsize=13, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend(fontsize=10)
ax1.set_ylim(0, 1.05)

# Plot 2: Security Level vs Sample Count
ax2.plot(m_range, security_bits, 'r-', linewidth=2)
ax2.axvline(x=optimal_samples[3], color='g', linestyle='--', linewidth=1.5,
            label=f'Optimal m={optimal_samples[3]}')
ax2.fill_between(m_range, 0, security_bits, alpha=0.3, color='red')
ax2.set_xlabel('Number of Samples (m)', fontsize=12)
ax2.set_ylabel('Security Level (bits)', fontsize=12)
ax2.set_title('Remaining Security Level vs Sample Count', fontsize=13, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend(fontsize=10)

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

print("2D plots generated successfully!")
print()

# Create 3D visualization
# Vary both dimension and sample count
dimensions = np.arange(30, 80, 5)
samples_3d = np.arange(40, 150, 5)
D, M = np.meshgrid(dimensions, samples_3d)

# Calculate success probability for each combination
Z = np.zeros_like(D, dtype=float)
for i in range(D.shape[0]):
    for j in range(D.shape[1]):
        n_val = D[i, j]
        m_val = M[i, j]
        Z[i, j] = attack_success_probability(n_val, q, sigma, m_val)

# Create 3D surface plot
fig = plt.figure(figsize=(14, 10))

# 3D Surface plot
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(D, M, Z, cmap='viridis', alpha=0.9, 
                        edgecolor='none', antialiased=True)
ax1.set_xlabel('Dimension (n)', fontsize=11, labelpad=10)
ax1.set_ylabel('Sample Count (m)', fontsize=11, labelpad=10)
ax1.set_zlabel('Success Probability', fontsize=11, labelpad=10)
ax1.set_title('3D Surface: Attack Success Probability', fontsize=12, fontweight='bold', pad=20)
ax1.view_init(elev=25, azim=45)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=5)

# 3D Contour plot
ax2 = fig.add_subplot(122, projection='3d')
contours = ax2.contour3D(D, M, Z, 50, cmap='plasma', alpha=0.8)
ax2.set_xlabel('Dimension (n)', fontsize=11, labelpad=10)
ax2.set_ylabel('Sample Count (m)', fontsize=11, labelpad=10)
ax2.set_zlabel('Success Probability', fontsize=11, labelpad=10)
ax2.set_title('3D Contour: Attack Success Probability', fontsize=12, fontweight='bold', pad=20)
ax2.view_init(elev=25, azim=45)
fig.colorbar(contours, ax=ax2, shrink=0.5, aspect=5)

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

print("3D plots generated successfully!")
print()

# Trade-off analysis
print("=== Trade-off Analysis ===")
print("Sample Count vs Security Level for P_target = 0.95:")
print("-" * 60)

trade_off_samples = [optimal_samples[3] - 20, optimal_samples[3], optimal_samples[3] + 20]
for m in trade_off_samples:
    if m > n:
        prob = attack_success_probability(n, q, sigma, m)
        sec = estimate_security_bits(n, q, sigma, m)
        efficiency = prob / m  # success probability per sample
        print(f"m = {m:3d}: P_success = {prob:.4f}, Security ≈ {sec:.1f} bits, Efficiency = {efficiency:.6f}")

print()
print("=== Conclusion ===")
print(f"For target success probability of 0.95:")
print(f"  - Minimum samples needed: {optimal_samples[3]}")
print(f"  - This is {optimal_samples[3] - n} samples above the theoretical minimum of n={n}")
print(f"  - Remaining security: ≈{estimate_security_bits(n, q, sigma, optimal_samples[3]):.1f} bits")
print()
print("Key insight: Adding more samples beyond the optimal point")
print("provides diminishing returns in success probability while")
print("continuously reducing the security of the system.")
print()
print("### Execution Results ###")
# Placeholder for execution results

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.linalg import qr
import time

def babai_algorithm(basis, target):
    """
    Babai's nearest plane algorithm for CVP approximation
    
    Args:
        basis: numpy array of shape (n, n) representing lattice basis
        target: numpy array of shape (n,) representing target point
    
    Returns:
        closest_point: approximation of closest lattice point
        coefficients: integer coefficients
    """
    # Gram-Schmidt orthogonalization
    Q, R = qr(basis.T)
    
    # Solve for coefficients in the orthogonal basis
    coeffs_real = np.linalg.solve(basis.T, target)
    
    # Round to nearest integers (Babai's rounding)
    coeffs_int = np.round(coeffs_real).astype(int)
    
    # Compute the closest lattice point
    closest_point = basis.T @ coeffs_int
    
    return closest_point, coeffs_int

def enumerate_cvp(basis, target, search_radius=5):
    """
    Brute force enumeration to find exact CVP solution
    
    Args:
        basis: numpy array of shape (n, n)
        target: numpy array of shape (n,)
        search_radius: how far to search in coefficient space
    
    Returns:
        best_point: exact closest lattice point
        best_coeffs: corresponding coefficients
        min_distance: minimum distance found
    """
    n = basis.shape[1]
    best_distance = float('inf')
    best_point = None
    best_coeffs = None
    
    # Generate all integer combinations within search radius
    ranges = [range(-search_radius, search_radius + 1) for _ in range(n)]
    
    import itertools
    for coeffs in itertools.product(*ranges):
        coeffs_array = np.array(coeffs)
        lattice_point = basis.T @ coeffs_array
        distance = np.linalg.norm(target - lattice_point)
        
        if distance < best_distance:
            best_distance = distance
            best_point = lattice_point
            best_coeffs = coeffs_array
    
    return best_point, best_coeffs, best_distance

def generate_lattice_points(basis, range_limit=5):
    """
    Generate lattice points for visualization
    """
    n = basis.shape[1]
    points = []
    
    import itertools
    ranges = [range(-range_limit, range_limit + 1) for _ in range(n)]
    
    for coeffs in itertools.product(*ranges):
        coeffs_array = np.array(coeffs)
        point = basis.T @ coeffs_array
        points.append(point)
    
    return np.array(points)

# Main execution
print("=" * 60)
print("Closest Vector Problem (CVP) Solver")
print("=" * 60)

# Define the 2D lattice basis
basis_2d = np.array([[3, 1],
                      [1, 2]])

target_2d = np.array([5.7, 3.2])

print("\nLattice Basis B:")
print(basis_2d)
print(f"\nTarget Point t: {target_2d}")

# Solve using Babai's algorithm
print("\n" + "=" * 60)
print("Method 1: Babai's Nearest Plane Algorithm (Fast Approximation)")
print("=" * 60)

start_time = time.time()
babai_point, babai_coeffs = babai_algorithm(basis_2d, target_2d)
babai_time = time.time() - start_time

babai_distance = np.linalg.norm(target_2d - babai_point)

print(f"Coefficients: {babai_coeffs}")
print(f"Closest Point (approx): {babai_point}")
print(f"Distance: {babai_distance:.6f}")
print(f"Computation Time: {babai_time:.6f} seconds")

# Solve using enumeration (exact solution)
print("\n" + "=" * 60)
print("Method 2: Exhaustive Enumeration (Exact Solution)")
print("=" * 60)

start_time = time.time()
exact_point, exact_coeffs, exact_distance = enumerate_cvp(basis_2d, target_2d, search_radius=10)
enum_time = time.time() - start_time

print(f"Coefficients: {exact_coeffs}")
print(f"Closest Point (exact): {exact_point}")
print(f"Distance: {exact_distance:.6f}")
print(f"Computation Time: {enum_time:.6f} seconds")

# 3D Example
print("\n" + "=" * 60)
print("3D Lattice Example")
print("=" * 60)

basis_3d = np.array([[4, 1, 0],
                      [1, 3, 1],
                      [0, 1, 3]])

target_3d = np.array([7.5, 5.2, 4.8])

print("\n3D Lattice Basis B:")
print(basis_3d)
print(f"\nTarget Point t: {target_3d}")

# Babai for 3D
print("\nBabai's Algorithm (3D):")
babai_point_3d, babai_coeffs_3d = babai_algorithm(basis_3d, target_3d)
babai_distance_3d = np.linalg.norm(target_3d - babai_point_3d)

print(f"Coefficients: {babai_coeffs_3d}")
print(f"Closest Point: {babai_point_3d}")
print(f"Distance: {babai_distance_3d:.6f}")

# Exact solution for 3D
print("\nExhaustive Enumeration (3D):")
start_time = time.time()
exact_point_3d, exact_coeffs_3d, exact_distance_3d = enumerate_cvp(basis_3d, target_3d, search_radius=5)
enum_time_3d = time.time() - start_time

print(f"Coefficients: {exact_coeffs_3d}")
print(f"Closest Point: {exact_point_3d}")
print(f"Distance: {exact_distance_3d:.6f}")
print(f"Computation Time: {enum_time_3d:.6f} seconds")

# Visualization
print("\n" + "=" * 60)
print("Generating Visualizations...")
print("=" * 60)

# Generate lattice points for 2D
lattice_points_2d = generate_lattice_points(basis_2d, range_limit=5)

# Generate lattice points for 3D
lattice_points_3d = generate_lattice_points(basis_3d, range_limit=3)

# Create comprehensive visualization
fig = plt.figure(figsize=(20, 12))

# 2D Plot - Overview
ax1 = fig.add_subplot(2, 3, 1)
ax1.scatter(lattice_points_2d[:, 0], lattice_points_2d[:, 1], 
           c='lightblue', s=50, alpha=0.6, label='Lattice Points')
ax1.scatter(target_2d[0], target_2d[1], 
           c='red', s=200, marker='*', label='Target', zorder=5)
ax1.scatter(exact_point[0], exact_point[1], 
           c='green', s=150, marker='s', label='Closest Point (Exact)', zorder=5)
ax1.scatter(babai_point[0], babai_point[1], 
           c='orange', s=150, marker='^', label='Babai Approximation', zorder=5)

# Draw basis vectors
origin = np.array([0, 0])
ax1.arrow(origin[0], origin[1], basis_2d[0, 0], basis_2d[1, 0], 
         head_width=0.3, head_length=0.3, fc='blue', ec='blue', linewidth=2, alpha=0.7)
ax1.arrow(origin[0], origin[1], basis_2d[0, 1], basis_2d[1, 1], 
         head_width=0.3, head_length=0.3, fc='purple', ec='purple', linewidth=2, alpha=0.7)

ax1.plot([target_2d[0], exact_point[0]], [target_2d[1], exact_point[1]], 
        'g--', linewidth=2, label=f'Distance: {exact_distance:.3f}')

ax1.set_xlabel('X', fontsize=12)
ax1.set_ylabel('Y', fontsize=12)
ax1.set_title('2D CVP: Lattice and Target Point', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.axis('equal')

# 2D Plot - Zoomed in
ax2 = fig.add_subplot(2, 3, 2)
nearby_points = lattice_points_2d[
    (np.abs(lattice_points_2d[:, 0] - target_2d[0]) < 5) & 
    (np.abs(lattice_points_2d[:, 1] - target_2d[1]) < 5)
]
ax2.scatter(nearby_points[:, 0], nearby_points[:, 1], 
           c='lightblue', s=100, alpha=0.8, label='Nearby Lattice Points')
ax2.scatter(target_2d[0], target_2d[1], 
           c='red', s=300, marker='*', label='Target', zorder=5)
ax2.scatter(exact_point[0], exact_point[1], 
           c='green', s=200, marker='s', label='Closest Point', zorder=5)

for point in nearby_points:
    distance = np.linalg.norm(point - target_2d)
    ax2.plot([target_2d[0], point[0]], [target_2d[1], point[1]], 
            'gray', alpha=0.3, linewidth=1)
    
ax2.plot([target_2d[0], exact_point[0]], [target_2d[1], exact_point[1]], 
        'g-', linewidth=3, label=f'Min Distance: {exact_distance:.3f}')

ax2.set_xlabel('X', fontsize=12)
ax2.set_ylabel('Y', fontsize=12)
ax2.set_title('2D CVP: Zoomed View', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.axis('equal')

# Distance comparison plot
ax3 = fig.add_subplot(2, 3, 3)
methods = ['Babai\n(Approx)', 'Enumeration\n(Exact)']
distances = [babai_distance, exact_distance]
colors = ['orange', 'green']
bars = ax3.bar(methods, distances, color=colors, alpha=0.7, edgecolor='black', linewidth=2)

for i, (bar, dist) in enumerate(zip(bars, distances)):
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height,
            f'{dist:.4f}',
            ha='center', va='bottom', fontsize=12, fontweight='bold')

ax3.set_ylabel('Distance to Target', fontsize=12)
ax3.set_title('2D CVP: Method Comparison', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3, axis='y')

# 3D Plot - Main view
ax4 = fig.add_subplot(2, 3, 4, projection='3d')
ax4.scatter(lattice_points_3d[:, 0], lattice_points_3d[:, 1], lattice_points_3d[:, 2],
           c='lightblue', s=30, alpha=0.4, label='Lattice Points')
ax4.scatter(target_3d[0], target_3d[1], target_3d[2],
           c='red', s=300, marker='*', label='Target', zorder=5)
ax4.scatter(exact_point_3d[0], exact_point_3d[1], exact_point_3d[2],
           c='green', s=200, marker='s', label='Closest Point', zorder=5)
ax4.plot([target_3d[0], exact_point_3d[0]], 
        [target_3d[1], exact_point_3d[1]], 
        [target_3d[2], exact_point_3d[2]],
        'g-', linewidth=3, label=f'Distance: {exact_distance_3d:.3f}')

# Draw basis vectors
origin_3d = np.array([0, 0, 0])
for i in range(3):
    ax4.quiver(origin_3d[0], origin_3d[1], origin_3d[2],
              basis_3d[0, i], basis_3d[1, i], basis_3d[2, i],
              arrow_length_ratio=0.1, linewidth=2, alpha=0.7)

ax4.set_xlabel('X', fontsize=12)
ax4.set_ylabel('Y', fontsize=12)
ax4.set_zlabel('Z', fontsize=12)
ax4.set_title('3D CVP: Lattice Structure', fontsize=14, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

# 3D Plot - Different angle
ax5 = fig.add_subplot(2, 3, 5, projection='3d')
nearby_points_3d = lattice_points_3d[
    (np.abs(lattice_points_3d[:, 0] - target_3d[0]) < 6) & 
    (np.abs(lattice_points_3d[:, 1] - target_3d[1]) < 6) &
    (np.abs(lattice_points_3d[:, 2] - target_3d[2]) < 6)
]
ax5.scatter(nearby_points_3d[:, 0], nearby_points_3d[:, 1], nearby_points_3d[:, 2],
           c='lightblue', s=60, alpha=0.6, label='Nearby Lattice Points')
ax5.scatter(target_3d[0], target_3d[1], target_3d[2],
           c='red', s=300, marker='*', label='Target', zorder=5)
ax5.scatter(exact_point_3d[0], exact_point_3d[1], exact_point_3d[2],
           c='green', s=200, marker='s', label='Closest Point', zorder=5)
ax5.plot([target_3d[0], exact_point_3d[0]], 
        [target_3d[1], exact_point_3d[1]], 
        [target_3d[2], exact_point_3d[2]],
        'g-', linewidth=3)

ax5.set_xlabel('X', fontsize=12)
ax5.set_ylabel('Y', fontsize=12)
ax5.set_zlabel('Z', fontsize=12)
ax5.set_title('3D CVP: Zoomed View', fontsize=14, fontweight='bold')
ax5.legend(fontsize=10)
ax5.view_init(elev=20, azim=45)
ax5.grid(True, alpha=0.3)

# 3D Distance comparison
ax6 = fig.add_subplot(2, 3, 6)
methods_3d = ['Babai\n(Approx)', 'Enumeration\n(Exact)']
distances_3d = [babai_distance_3d, exact_distance_3d]
colors_3d = ['orange', 'green']
bars_3d = ax6.bar(methods_3d, distances_3d, color=colors_3d, alpha=0.7, edgecolor='black', linewidth=2)

for i, (bar, dist) in enumerate(zip(bars_3d, distances_3d)):
    height = bar.get_height()
    ax6.text(bar.get_x() + bar.get_width()/2., height,
            f'{dist:.4f}',
            ha='center', va='bottom', fontsize=12, fontweight='bold')

ax6.set_ylabel('Distance to Target', fontsize=12)
ax6.set_title('3D CVP: Method Comparison', fontsize=14, fontweight='bold')
ax6.grid(True, alpha=0.3, axis='y')

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

print("\nVisualization complete!")
print("=" * 60)