fpylll: A Python Interface to the fplll Lattice Reduction Library

# Import fpylll and related libraries
import fpylll
from fpylll import IntegerMatrix, GSO, LLL, BKZ, FPLLL
from fpylll.algorithms.bkz import BKZReduction
from fpylll.tools.quality import basis_quality
import numpy as np
import matplotlib.pyplot as plt

# Verify installation
print(f"fpylll version: {fpylll.__version__}")
print(f"Available float types: {fpylll.config.float_types}")

# Set random seed for reproducibility
FPLLL.set_random_seed(42)
print("fpylll successfully imported and configured!")

# 1. Create empty matrix and populate manually
A = IntegerMatrix(3, 3)
A[0, 0] = 1; A[0, 1] = 2; A[0, 2] = 3
A[1, 0] = 4; A[1, 1] = 5; A[1, 2] = 6  
A[2, 0] = 7; A[2, 1] = 8; A[2, 2] = 9
print("Manual construction:")
print(A)

# 2. Create from Python lists/arrays
B = IntegerMatrix.from_matrix([[1, 0, 2], [0, 1, 3], [0, 0, 1]])
print("\nFrom Python list:")
print(B)

# 3. Create identity matrix
I = IntegerMatrix.identity(4)
print("\nIdentity matrix:")
print(I)

# 4. Matrix properties
print(f"\nMatrix dimensions: {A.nrows} × {A.ncols}")
print(f"Matrix norm of first row: {A[0].norm()}")
print(f"Matrix determinant estimate: {A.get_max_exp()}")

# Generate different types of random lattices
print("=== Random Lattice Generation Examples ===\n")

# 1. Uniform random matrix
uniform_lattice = IntegerMatrix.random(4, "uniform", bits=10)
print("1. Uniform random lattice (4×4, 10-bit entries):")
print(uniform_lattice)

# 2. NTRU-like lattice  
ntru_lattice = IntegerMatrix.random(5, "ntrulike", q=127)
print("\n2. NTRU-like lattice (5×10, q=127):")
print(ntru_lattice)

# 3. q-ary lattice (SIS/LWE related)
qary_lattice = IntegerMatrix.random(6, "qary", q=127, k=3)
print("\n3. q-ary lattice (6×6, q=127, k=3):")
print(qary_lattice)

# 4. Integer relation lattice
intrel_lattice = IntegerMatrix.random(4, "intrel", bits=20)
print("\n4. Integer relation lattice (4×5, 20-bit):")
print(intrel_lattice)

# 5. Compute and display basic properties
def analyze_lattice(A, name):
    print(f"\n--- Analysis of {name} ---")
    print(f"Dimensions: {A.nrows}×{A.ncols}")
    print(f"First vector norm: {A[0].norm():.2f}")
    print(f"Last vector norm: {A[A.nrows-1].norm():.2f}")
    if A.nrows == A.ncols:
        # Estimate volume for square matrices
        volume_est = 1
        for i in range(A.nrows):
            volume_est *= A[i].norm()
        print(f"Volume estimate: {volume_est**(1/A.nrows):.2f}")

analyze_lattice(uniform_lattice, "Uniform Lattice")
analyze_lattice(ntru_lattice, "NTRU Lattice")

# Create a lattice basis and GSO object
A = IntegerMatrix.random(5, "uniform", bits=10)
print("Original basis:")
print(A)

# Create GSO object and compute orthogonalization
M = GSO.Mat(A)
M.update_gso()

print("\n=== Gram-Schmidt Analysis ===")

# Access μ coefficients  
print("\nμ coefficients (lower triangular):")
for i in range(1, A.nrows):
    for j in range(i):
        mu_ij = M.get_mu(i, j)
        print(f"μ_{i},{j} = {mu_ij:.4f}", end="  ")
    print()

# Access r coefficients (squared norms of b_i*)
print("\nr coefficients (squared norms):")
for i in range(A.nrows):
    r_ii = M.get_r(i, i)
    print(f"r_{i},{i} = ||b*_{i}||² = {r_ii:.2f}")

# Gram matrix elements
print("\nGram matrix (first 3×3):")
for i in range(min(3, A.nrows)):
    for j in range(min(3, A.ncols)):
        if i >= j:
            gram_ij = M.get_gram(i, j)
            print(f"{gram_ij:8.0f}", end=" ")
        else:
            print(f"{'':8s}", end=" ")
    print()

# Quality measures
print(f"\nLog determinant: {M.get_log_det(0, A.nrows):.4f}")
print(f"Root determinant: {M.get_root_det(0, A.nrows):.4f}")

# Create a random lattice for LLL demonstration
print("=== LLL Reduction Example ===\n")

# Generate a challenging lattice (larger entries, less structured)
A = IntegerMatrix.random(8, "uniform", bits=30)
print("Original basis vectors norms:")
for i in range(A.nrows):
    print(f"||b_{i}|| = {A[i].norm():.2f}")

# Compute initial quality measures
def compute_quality(B, name=""):
    """Compute and display basis quality measures"""
    M = GSO.Mat(B)
    M.update_gso()
    
    # Compute Hermite factor
    det = M.get_root_det(0, B.nrows)
    hermite_factor = B[0].norm() / (det ** (1.0/B.nrows))
    
    # Compute orthogonality defect
    volume_basis = 1.0
    volume_gso = 1.0
    for i in range(B.nrows):
        volume_basis *= B[i].norm()
        volume_gso *= M.get_r(i, i) ** 0.5
    
    orthogonality_defect = volume_basis / volume_gso
    
    print(f"\n{name} Quality Measures:")
    print(f"  First vector norm: {B[0].norm():.2f}")
    print(f"  Root determinant: {det:.2f}")
    print(f"  Hermite factor: {hermite_factor:.4f}")
    print(f"  Orthogonality defect: {orthogonality_defect:.2f}")
    
    return hermite_factor, orthogonality_defect

# Measure quality before reduction
hf_before, od_before = compute_quality(A, "BEFORE LLL")

# Perform LLL reduction
print(f"\nPerforming LLL reduction with δ=0.99, η=0.51...")
A_copy = IntegerMatrix(A)  # Make a copy
LLL.reduction(A_copy)

# Measure quality after reduction  
hf_after, od_after = compute_quality(A_copy, "AFTER LLL")

# Show improvement
print(f"\n=== Improvement ===")
print(f"Hermite factor: {hf_before:.4f} → {hf_after:.4f} (factor {hf_before/hf_after:.2f})")
print(f"Orthogonality defect: {od_before:.2f} → {od_after:.2f} (factor {od_before/od_after:.2f})")

# Verify LLL-reducedness
is_reduced = LLL.is_reduced(A_copy)
print(f"Is LLL-reduced: {is_reduced}")

print("\nReduced basis vector norms:")
for i in range(A_copy.nrows):
    print(f"||b_{i}|| = {A_copy[i].norm():.2f}")

# BKZ reduction example with different block sizes
print("=== BKZ Reduction Comparison ===\n")

# Create a more challenging lattice 
A = IntegerMatrix.random(20, "qary", q=1024, k=10)
print(f"Generated {A.nrows}×{A.ncols} q-ary lattice with q=1024, k=10")

# Helper function for BKZ analysis
def run_bkz_analysis(basis, block_size, max_tours=10):
    """Run BKZ and analyze results"""
    from fpylll.algorithms.bkz2 import BKZReduction as BKZ2
    
    # Create copy for reduction
    B = IntegerMatrix(basis)
    
    # First apply LLL for fair comparison
    LLL.reduction(B)
    
    print(f"\n--- BKZ-{block_size} Analysis ---")
    
    # Measure initial quality after LLL
    M_before = GSO.Mat(B)
    M_before.update_gso()
    det = M_before.get_root_det(0, B.nrows)
    hf_before = B[0].norm() / (det ** (1.0/B.nrows))
    
    print(f"After LLL preprocessing:")
    print(f"  First vector norm: {B[0].norm():.2f}")
    print(f"  Hermite factor: {hf_before:.6f}")
    
    if block_size > 2:
        # Apply BKZ reduction
        from fpylll import BKZ as BKZ_param
        param = BKZ_param.Param(
            block_size=block_size,
            max_loops=max_tours,
            strategies=BKZ_param.DEFAULT_STRATEGY
        )
        
        bkz = BKZ2(GSO.Mat(B))
        bkz(param)
    
    # Measure final quality
    M_after = GSO.Mat(B) 
    M_after.update_gso()
    hf_after = B[0].norm() / (det ** (1.0/B.nrows))
    
    print(f"After BKZ-{block_size}:")
    print(f"  First vector norm: {B[0].norm():.2f}")
    print(f"  Hermite factor: {hf_after:.6f}")
    print(f"  Improvement factor: {hf_before/hf_after:.3f}")
    
    return hf_after, B

# Compare different block sizes
block_sizes = [2, 10, 20]  # Start with smaller block sizes
results = {}

for beta in block_sizes:
    try:
        hf, reduced_basis = run_bkz_analysis(A, beta)
        results[beta] = hf
    except Exception as e:
        print(f"BKZ-{beta} failed: {e}")
        results[beta] = float('inf')

print(f"\n=== BKZ Comparison Summary ===")
print("Block Size | Hermite Factor")
print("-" * 25)
for beta, hf in results.items():
    if hf != float('inf'):
        print(f"BKZ-{beta:2d}     | {hf:.6f}")
    else:
        print(f"BKZ-{beta:2d}     | Failed")

# Theoretical Hermite factors for comparison
print(f"\nTheoretical asymptotic Hermite factors:")
for beta in block_sizes:
    if beta == 2:
        theoretical_hf = float(2**(1/4))  # LLL bound
    else:
        # Simplified BKZ bound - convert to float for formatting
        import math
        theoretical_hf = float((beta/(2*math.pi*math.e) * (math.pi*beta)**(1/beta))**(1/(2*(beta-1))))
    print(f"BKZ-{beta}: δ ≈ {theoretical_hf:.6f}")

# SVP and CVP examples
print("=== SVP and CVP Examples ===\n")

# Create a small lattice for exact SVP solving
A = IntegerMatrix.random(6, "uniform", bits=8)
LLL.reduction(A)  # Preprocess with LLL

print("Working with 6-dimensional lattice:")
print(A)

# 1. Babai's algorithm for approximate CVP
print(f"\n1. Babai's Nearest Plane Algorithm")
M = GSO.Mat(A)
M.update_gso()

# Choose a random target vector
target = [100, 50, 75, 25, 125, 80]
print(f"Target vector: {target}")

# Apply Babai's algorithm
closest_coeff = M.babai(target)
closest_vector = A.multiply_left(closest_coeff)

print(f"Closest lattice vector: {closest_vector}")
print(f"Distance: {sum((target[i] - closest_vector[i])**2 for i in range(len(target)))**0.5:.2f}")

# 2. Using built-in CVP solver
from fpylll import CVP
closest_vector_cvp = CVP.babai(A, target)
print(f"CVP.babai result: {closest_vector_cvp}")

# 3. Enumeration for small lattices
print(f"\n2. Enumeration Example")
from fpylll import Enumeration

# Create enumeration object
enum = Enumeration(M)

# Set enumeration radius (use first vector as bound)
radius = A[0].norm()**2
print(f"Enumeration radius: {radius**.5:.2f}")

try:
    # Find short vectors (this may be computationally intensive)
    # We use a smaller radius to find multiple solutions
    solutions = enum.enumerate(0, A.nrows, 0.8 * radius, 0)
    
    print(f"Found {len(solutions)} vectors within radius {(0.8 * radius)**0.5:.2f}")
    
    # Display shortest vectors found
    if solutions:
        solutions.sort(key=lambda x: x[1])  # Sort by length
        print("Shortest vectors found:")
        for i, (coeff, length) in enumerate(solutions[:3]):
            vector = A.multiply_left(coeff)
            print(f"  Vector {i+1}: {vector}, length: {length**0.5:.2f}")
            
except Exception as e:
    print(f"Enumeration failed or took too long: {e}")

# 4. Gaussian Heuristic prediction
print(f"\n3. Gaussian Heuristic Analysis")
from fpylll.util import gaussian_heuristic

# Compute Gaussian heuristic prediction
gh_length = gaussian_heuristic(M.r())**0.5
print(f"Gaussian heuristic prediction: {gh_length:.2f}")
print(f"Actual first vector length: {A[0].norm():.2f}")
print(f"Ratio (actual/predicted): {A[0].norm()/gh_length:.3f}")

# 5. Volume and determinant analysis
det = M.get_root_det(0, A.nrows)
print(f"Lattice determinant: {det:.2e}")
print(f"Volume per lattice point: {det:.2e}")

# Breaking a subset sum (knapsack) problem using lattice reduction

# Generate a subset sum problem instance
import random
random.seed(int(42))  # Explicitly convert to Python int

# Parameters for the subset sum problem
n = int(20)  # Explicitly convert to Python int
density = float(0.9)  # Density of the knapsack (affects difficulty)

# Generate random weights
weights = [random.randint(1, 2**(int(n*density))) for _ in range(n)]
print(f"Weights: {weights[:5]}...")  # Show first 5 weights

# Generate a random subset and compute the target sum
subset_indices = random.sample(range(n), n//2)
target = sum(weights[i] for i in subset_indices)
print(f"Target sum: {target}")
print(f"Subset indices: {subset_indices}")

# Construct the lattice for subset sum attack
# We use the standard lattice construction for subset sum
from fpylll import IntegerMatrix

# Create the lattice matrix
# The lattice has dimension (n+1) x (n+1)
M = IntegerMatrix(n+1, n+1)

# Fill the matrix according to the subset sum lattice construction
scaling_factor = int(2)  # Scaling factor to balance the lattice

for i in range(n):
    M[i, i] = int(2)  # Identity part scaled
    M[i, n] = int(weights[i])  # Weights in the last column

M[n, n] = int(target)  # Target in the bottom-right corner

print(f"Lattice dimension: {M.nrows} x {M.ncols}")
# Note: determinant computation omitted for efficiency

# Apply BKZ reduction to find short vectors
from fpylll import BKZ

# Use a moderate block size for efficiency
block_size = int(10)
param = BKZ.Param(block_size=block_size, strategies=BKZ.DEFAULT_STRATEGY)

# Reduce the lattice
print(f"\nApplying BKZ-{block_size} reduction...")
BKZ.reduction(M, param)

# Check if we found a solution
print("\nChecking for solutions in reduced basis:")
for i in range(min(int(5), M.nrows)):  # Check first few vectors
    vector = [M[i, j] for j in range(M.ncols)]
    
    # Check if this vector corresponds to a subset sum solution
    if abs(vector[-1]) <= int(1):  # The last coordinate should be 0 or ±1
        # Extract the subset from the vector
        subset_vector = vector[:-1]
        
        # Check if this gives a valid subset (binary vector after scaling)
        if all(abs(x) <= int(1) for x in subset_vector):
            # Compute the corresponding subset sum
            computed_sum = sum(int(weights[j]) * abs(subset_vector[j]) for j in range(n))
            
            if abs(computed_sum - target) <= abs(vector[-1]):
                print(f"Found solution vector: {vector}")
                solution_indices = [j for j in range(n) if abs(subset_vector[j]) > int(0)]
                print(f"Solution subset indices: {solution_indices}")
                print(f"Original subset indices: {subset_indices}")
                print(f"Match: {set(solution_indices) == set(subset_indices)}")
                break
else:
    print("No clear solution found in the first few vectors")

# Demonstration of fpylll utility functions

from fpylll import IntegerMatrix, LLL, BKZ, GSO
from fpylll.util import gaussian_heuristic
import math

# Generate different types of random lattices
print("=== Random Lattice Generation ===")

# 1. Uniform random lattice
n = 30
A = IntegerMatrix.random(n, "uniform", bits=30)
print(f"Uniform random lattice: {A.nrows}x{A.ncols}, bits=30")

# 2. q-ary lattice (important in cryptography)
q = 97  # Small prime
m = 40  # Number of samples
A_qary = IntegerMatrix.random(m, "qary", k=n, bits=30, q=q)
print(f"q-ary lattice: q={q}, dimension={A_qary.nrows}x{A_qary.ncols}")

# 3. NTRU-like lattice
A_ntru = IntegerMatrix.random(n, "ntrulike", bits=10, q=q)
print(f"NTRU-like lattice: {A_ntru.nrows}x{A_ntru.ncols}, q={q}")

print("\n=== Quality Measures ===")

# Work with a small example for clear demonstration
B = IntegerMatrix.random(15, "uniform", bits=20)
print(f"Original lattice dimension: {B.nrows}x{B.ncols}")

# Apply LLL reduction
LLL.reduction(B)

# Create GSO object to compute quality measures
M = GSO.Mat(B)
M.update_gso()

# Compute various quality measures
n = B.nrows
vol = B.det()

# First vector length
first_length = sum(B[0, j]**2 for j in range(B.ncols))**0.5
print(f"First vector length: {first_length:.3f}")

# Theoretical shortest vector (Gaussian heuristic)
gh = gaussian_heuristic([M.get_r(i, i) for i in range(n)])
print(f"Gaussian heuristic: {gh:.3f}")

# Hermite factor
hermite_factor = (first_length / abs(vol)**(1.0/n))**n
print(f"Hermite factor: {hermite_factor:.6f}")

# Root Hermite factor  
root_hermite = hermite_factor**(1.0/n)
print(f"Root Hermite factor: {root_hermite:.6f}")

# Orthogonality defect
product_lengths = 1.0
for i in range(n):
    length = sum(B[i, j]**2 for j in range(B.ncols))**0.5
    product_lengths *= length

orthogonality_defect = product_lengths / abs(vol)
print(f"Orthogonality defect: {orthogonality_defect:.3e}")

print("\n=== Precision Options ===")

# Demonstrate different precision options
from fpylll import FPLLL

# Check current precision
print(f"Current default precision: {FPLLL.get_precision()} bits")

# You can set higher precision for more accurate computations
# FPLLL.set_precision(100)  # Uncomment to set 100-bit precision
# print(f"New precision: {FPLLL.get_precision()} bits")

# Different precision modes affect GSO computations
print("GSO with different floating point types:")
print("- FP_NR<double>: Fast, double precision")
print("- FP_NR<long double>: Extended precision") 
print("- FP_NR<mpfr_t>: Arbitrary precision")

print("\n=== Performance Comparison ===")

import time

# Compare LLL vs BKZ performance on same lattice
test_lattice = IntegerMatrix.random(25, "uniform", bits=25)

# LLL timing
test_copy = IntegerMatrix.from_matrix(test_lattice)
start_time = time.time()
LLL.reduction(test_copy)
lll_time = time.time() - start_time

# BKZ-10 timing  
test_copy = IntegerMatrix.from_matrix(test_lattice)
start_time = time.time()
BKZ.reduction(test_copy, BKZ.Param(block_size=10))
bkz10_time = time.time() - start_time

# BKZ-20 timing
test_copy = IntegerMatrix.from_matrix(test_lattice)
start_time = time.time()
BKZ.reduction(test_copy, BKZ.Param(block_size=20))
bkz20_time = time.time() - start_time

print(f"LLL reduction time: {lll_time:.4f} seconds")
print(f"BKZ-10 reduction time: {bkz10_time:.4f} seconds")  
print(f"BKZ-20 reduction time: {bkz20_time:.4f} seconds")
print(f"BKZ-20 vs LLL slowdown: {bkz20_time/lll_time:.1f}x")

# Example: Coppersmith's Method for Small Roots

# This demonstrates the basic idea behind Coppersmith's method
# We'll solve f(x) = 0 mod N where we know x is small

from fpylll import IntegerMatrix, LLL
import math

def coppersmith_univariate_example():
    """
    Example of using lattice reduction to find small roots of 
    a polynomial modulo a large integer (simplified Coppersmith)
    """
    
    # Example: f(x) = x^2 + a*x + b mod N
    # We want to find small root x₀ such that f(x₀) ≡ 0 (mod N)
    
    # Parameters
    N = 1009 * 1013  # Small composite for demonstration
    a = 123
    b = 456
    
    # The polynomial f(x) = x^2 + 123*x + 456
    # Let's say we know the root is small (< 100)
    
    print(f"Finding small roots of f(x) = x² + {a}x + {b} mod {N}")
    
    # Generate a small root for testing
    true_root = 37  # This should be unknown in practice
    check = (true_root**2 + a*true_root + b) % N
    print(f"Verification: f({true_root}) ≡ {check} (mod {N})")
    
    # Coppersmith's method construction
    # We construct a lattice where short vectors correspond to
    # polynomials that have the same small root but are zero over integers
    
    d = 2  # Degree of polynomial
    m = 2  # Parameter controlling lattice dimension
    X = 100  # Bound on the root size
    
    # Construct lattice for Coppersmith's method
    # This is a simplified version - real implementation needs more care
    
    dim = m + d  # Lattice dimension
    L = IntegerMatrix(dim, dim)
    
    # Fill the lattice matrix
    # Each row corresponds to a polynomial of the form x^i * f(x)^j * N^(m-j)
    row = 0
    
    for j in range(m):
        for i in range(d):
            if row < dim:
                # Polynomial: x^i * f(x)^j * N^(m-j)
                # Coefficient of x^k in this polynomial
                L[row, i] = X**i * N**(m-j)  # Simplified
                row += 1
    
    # Add some structure to make it more realistic
    for i in range(dim):
        if L[i, i] == 0:
            L[i, i] = X**i
    
    print(f"Lattice dimension: {dim}x{dim}")
    print(f"Original determinant: {abs(L.det())}")
    
    # Apply LLL reduction
    LLL.reduction(L)
    print(f"Reduced determinant: {abs(L.det())}")
    
    # In a real implementation, we would:
    # 1. Extract polynomials from short lattice vectors
    # 2. Compute their GCD over integers
    # 3. Factor the resulting polynomial to find roots
    
    print("In practice, short vectors give polynomials with small integer roots")
    print("These can be factored to recover the original small root")
    
    return true_root

def simple_rsa_example():
    """
    Simplified example showing how lattices can attack RSA with partial information
    """
    print("\n=== RSA Partial Key Exposure Example ===")
    
    # Very small RSA parameters for demonstration
    p, q = 61, 67
    N = p * q
    phi = (p-1) * (q-1)
    e = 17  # Public exponent
    
    # Compute private exponent
    d = pow(e, -1, phi)  # Modular inverse
    
    print(f"RSA parameters: N={N}, e={e}, d={d}")
    
    # Suppose we know some bits of d (partial key exposure)
    # In practice, this might come from side-channel attacks
    
    known_bits = 6  # We know the lower 6 bits of d
    d_low = d & ((1 << known_bits) - 1)  # Extract lower bits
    d_high = d >> known_bits  # Unknown high bits
    
    print(f"Known low bits of d: {d_low} (correct: {d & ((1 << known_bits) - 1)})")
    print(f"Unknown high bits: {d_high}")
    
    # Set up lattice to find the unknown high bits
    # This uses the relation: e*d ≡ 1 (mod φ(N))
    # So e*d = 1 + k*φ(N) for some integer k
    
    # In practice, we'd use Coppersmith's method here
    # For now, we just verify our approach would work
    
    # The equation becomes: e*(d_high * 2^known_bits + d_low) = 1 + k*φ(N)
    # Rearranging: e*d_high * 2^known_bits = 1 + k*φ(N) - e*d_low
    
    scale = 2**known_bits
    target = 1 - e * d_low  # This should be ≡ k*φ(N) (mod something)
    
    print(f"Attack would solve for d_high such that:")
    print(f"e * d_high * {scale} ≡ {target} + k*φ(N)")
    print(f"This can be solved using Coppersmith's method with appropriate lattice")
    
    return d

# Run the examples
print("=== Coppersmith Method Demonstration ===")
root = coppersmith_univariate_example()
private_key = simple_rsa_example()

print(f"\n=== Summary ===")
print("These examples show the basic principles behind lattice-based cryptanalysis:")
print("1. Transform the problem into finding small solutions")
print("2. Construct a lattice where short vectors encode the solutions") 
print("3. Use LLL/BKZ to find short vectors")
print("4. Extract and verify the solutions")

# Performance optimization examples

from fpylll import IntegerMatrix, LLL, BKZ, GSO
import time

def benchmark_algorithms():
    """Compare performance of different reduction algorithms"""
    
    print("=== Performance Benchmarking ===")
    
    # Test on various lattice dimensions
    dimensions = [20, 30, 40, 50]
    
    for n in dimensions:
        print(f"\nDimension {n}:")
        
        # Generate test lattice
        A = IntegerMatrix.random(n, "uniform", bits=30)
        
        # LLL benchmark
        A_copy = IntegerMatrix.from_matrix(A)
        start = time.time()
        LLL.reduction(A_copy)
        lll_time = time.time() - start
        
        # Compute quality after LLL
        M_lll = GSO.Mat(A_copy)
        M_lll.update_gso()
        lll_quality = M_lll.get_r(0, 0)**0.5 / (abs(A_copy.det())**(1.0/n))
        
        # BKZ-10 benchmark
        A_copy = IntegerMatrix.from_matrix(A)
        start = time.time()
        BKZ.reduction(A_copy, BKZ.Param(block_size=10))
        bkz10_time = time.time() - start
        
        # Compute quality after BKZ-10
        M_bkz = GSO.Mat(A_copy)
        M_bkz.update_gso()
        bkz10_quality = M_bkz.get_r(0, 0)**0.5 / (abs(A_copy.det())**(1.0/n))
        
        # BKZ-20 benchmark (only for smaller dimensions)
        if n <= 40:
            A_copy = IntegerMatrix.from_matrix(A)
            start = time.time()
            BKZ.reduction(A_copy, BKZ.Param(block_size=20))
            bkz20_time = time.time() - start
            
            M_bkz20 = GSO.Mat(A_copy)
            M_bkz20.update_gso()
            bkz20_quality = M_bkz20.get_r(0, 0)**0.5 / (abs(A_copy.det())**(1.0/n))
            
            print(f"  LLL:    {lll_time:.3f}s, quality: {lll_quality:.4f}")
            print(f"  BKZ-10: {bkz10_time:.3f}s, quality: {bkz10_quality:.4f}")
            print(f"  BKZ-20: {bkz20_time:.3f}s, quality: {bkz20_quality:.4f}")
        else:
            print(f"  LLL:    {lll_time:.3f}s, quality: {lll_quality:.4f}")
            print(f"  BKZ-10: {bkz10_time:.3f}s, quality: {bkz10_quality:.4f}")
            print(f"  BKZ-20: Skipped (too slow for demo)")

def demonstrate_strategies():
    """Show different BKZ strategies and their effects"""
    
    print("\n=== BKZ Strategy Comparison ===")
    
    # Create test lattice
    n = 30
    A = IntegerMatrix.random(n, "uniform", bits=25)
    
    # Default strategy
    A_copy = IntegerMatrix.from_matrix(A)
    start = time.time()
    param_default = BKZ.Param(block_size=15, strategies=BKZ.DEFAULT_STRATEGY)
    BKZ.reduction(A_copy, param_default)
    default_time = time.time() - start
    
    M_default = GSO.Mat(A_copy)
    M_default.update_gso()
    default_quality = M_default.get_r(0, 0)**0.5
    
    print(f"Default strategy: {default_time:.3f}s, ||b₁|| = {default_quality:.2f}")
    
    # Custom strategy with more enumeration
    A_copy = IntegerMatrix.from_matrix(A)
    start = time.time()
    
    # Create custom strategy (simplified example)
    param_custom = BKZ.Param(block_size=15)
    param_custom.strategies = BKZ.DEFAULT_STRATEGY  # In practice, you'd customize this
    
    BKZ.reduction(A_copy, param_custom)
    custom_time = time.time() - start
    
    M_custom = GSO.Mat(A_copy)
    M_custom.update_gso()
    custom_quality = M_custom.get_r(0, 0)**0.5
    
    print(f"Custom strategy:  {custom_time:.3f}s, ||b₁|| = {custom_quality:.2f}")

def precision_examples():
    """Demonstrate different precision settings"""
    
    print("\n=== Precision Settings ===")
    
    from fpylll import FPLLL
    
    # Show current precision
    current_prec = FPLLL.get_precision()
    print(f"Current precision: {current_prec} bits")
    
    # For very ill-conditioned lattices, higher precision helps
    # Create an ill-conditioned example
    n = 15
    A = IntegerMatrix(n, n)
    
    # Create identity matrix with exponentially growing diagonal
    for i in range(n):
        A[i, i] = 2**(i*2)  # Exponentially growing
    
    print(f"Condition number estimate: ~2^{(n-1)*2}")
    
    # Standard precision reduction
    A_copy = IntegerMatrix.from_matrix(A)
    LLL.reduction(A_copy)
    
    M = GSO.Mat(A_copy)
    M.update_gso()
    standard_quality = M.get_r(0, 0)**0.5
    
    print(f"Standard precision result: ||b₁|| = {standard_quality:.3e}")
    
    # Note: In practice, you might set higher precision for ill-conditioned cases
    # FPLLL.set_precision(100)  # Uncomment for higher precision
    # Then repeat the reduction...
    
    print("For ill-conditioned lattices, consider:")
    print("- Using higher precision (FPLLL.set_precision)")
    print("- Checking numerical stability")
    print("- Using exact arithmetic when possible")

def memory_optimization():
    """Tips for memory-efficient lattice reduction"""
    
    print("\n=== Memory Optimization Tips ===")
    
    # Memory usage scales quadratically with dimension for GSO
    n = 50
    A = IntegerMatrix.random(n, "uniform", bits=20)
    
    print(f"For dimension {n}:")
    print(f"- GSO matrix requires ~{n*n*8} bytes for double precision")
    print(f"- IntegerMatrix requires ~{n*n*8} bytes for the basis")
    print(f"- BKZ enumeration may need additional temporary storage")
    
    # Memory-efficient approach: avoid keeping multiple GSO objects
    print("\nMemory-efficient practices:")
    print("1. Reuse GSO objects when possible")
    print("2. Use in-place operations")
    print("3. Consider block-wise processing for very large lattices")
    print("4. Free intermediate objects explicitly if needed")
    
    # Example of reusing GSO object
    M = GSO.Mat(A)
    M.update_gso()
    
    # Work with the GSO object...
    quality_before = M.get_r(0, 0)**0.5
    
    # Apply reduction (this modifies A in-place)
    BKZ.reduction(A, BKZ.Param(block_size=10))
    
    # Update the same GSO object
    M.update_gso()
    quality_after = M.get_r(0, 0)**0.5
    
    print(f"\nReusing GSO object:")
    print(f"Quality before BKZ: {quality_before:.3f}")
    print(f"Quality after BKZ:  {quality_after:.3f}")

# Run all optimization examples
benchmark_algorithms()
demonstrate_strategies()
precision_examples()
memory_optimization()

# Final comprehensive example: Complete lattice analysis workflow

from fpylll import IntegerMatrix, LLL, BKZ, GSO, Enumeration
from fpylll.util import gaussian_heuristic
import time
import math

def complete_lattice_analysis(dimension=25, bits=30):
    """
    Complete workflow demonstrating all major fpylll capabilities
    """
    print("=" * 60)
    print(f"COMPLETE LATTICE ANALYSIS WORKFLOW")
    print("=" * 60)
    
    # Step 1: Generate and analyze initial lattice
    print(f"\n1. LATTICE GENERATION (dimension {dimension}, {bits} bits)")
    print("-" * 40)
    
    A = IntegerMatrix.random(dimension, "uniform", bits=bits)
    print(f"Generated {A.nrows}×{A.ncols} random lattice")
    print(f"Determinant: {A.det()}")
    
    # Initial quality analysis
    M = GSO.Mat(A)
    M.update_gso()
    
    initial_length = M.get_r(0, 0)**0.5
    vol = abs(A.det())
    initial_hermite = (initial_length / vol**(1.0/dimension))**dimension
    
    print(f"Initial first vector length: {initial_length:.3f}")
    print(f"Initial Hermite factor: {initial_hermite:.6f}")
    
    # Step 2: LLL Reduction
    print(f"\n2. LLL REDUCTION")
    print("-" * 40)
    
    start_time = time.time()
    LLL.reduction(A)
    lll_time = time.time() - start_time
    
    M.update_gso()
    lll_length = M.get_r(0, 0)**0.5
    lll_hermite = (lll_length / vol**(1.0/dimension))**dimension
    
    print(f"LLL time: {lll_time:.4f} seconds")
    print(f"LLL first vector length: {lll_length:.3f}")
    print(f"LLL Hermite factor: {lll_hermite:.6f}")
    print(f"Improvement factor: {initial_length/lll_length:.2f}")
    
    # Step 3: BKZ Reduction
    print(f"\n3. BKZ REDUCTION")
    print("-" * 40)
    
    block_size = min(20, dimension//2)
    print(f"Applying BKZ-{block_size}...")
    
    start_time = time.time()
    param = BKZ.Param(block_size=block_size, strategies=BKZ.DEFAULT_STRATEGY)
    BKZ.reduction(A, param)
    bkz_time = time.time() - start_time
    
    M.update_gso()
    bkz_length = M.get_r(0, 0)**0.5
    bkz_hermite = (bkz_length / vol**(1.0/dimension))**dimension
    
    print(f"BKZ-{block_size} time: {bkz_time:.4f} seconds")
    print(f"BKZ first vector length: {bkz_length:.3f}")
    print(f"BKZ Hermite factor: {bkz_hermite:.6f}")
    print(f"Improvement over LLL: {lll_length/bkz_length:.2f}")
    
    # Step 4: Quality Analysis
    print(f"\n4. QUALITY ANALYSIS")
    print("-" * 40)
    
    # Gaussian heuristic
    gh = gaussian_heuristic([M.get_r(i, i) for i in range(dimension)])
    print(f"Gaussian heuristic: {gh:.3f}")
    print(f"First vector vs GH: {bkz_length/gh:.3f}")
    
    # Orthogonality defect
    product_lengths = 1.0
    for i in range(dimension):
        row_length = sum(A[i, j]**2 for j in range(A.ncols))**0.5
        product_lengths *= row_length
    
    orthogonality_defect = product_lengths / vol
    print(f"Orthogonality defect: {orthogonality_defect:.3e}")
    
    # Successive minima approximation
    print(f"Successive minima approximation:")
    for i in range(min(5, dimension)):
        length_i = M.get_r(i, i)**0.5
        print(f"  λ_{i+1} ≈ {length_i:.3f}")
    
    # Step 5: Short Vector Search (for small dimensions)
    if dimension <= 20:
        print(f"\n5. ENUMERATION (SVP)")
        print("-" * 40)
        
        try:
            # Enumerate short vectors
            start_time = time.time()
            enum = Enumeration(M)
            
            # Find vectors within 1.1 times the first vector length
            radius = bkz_length * 1.1
            solutions = enum.enumerate(0, dimension, radius**2, 0)
            enum_time = time.time() - start_time
            
            print(f"Enumeration time: {enum_time:.4f} seconds")
            print(f"Search radius: {radius:.3f}")
            print(f"Vectors found: {len(solutions)}")
            
            if solutions:
                shortest = min(solutions, key=lambda x: x[0])
                shortest_length = shortest[0]**0.5
                print(f"Shortest vector length: {shortest_length:.3f}")
                print(f"Ratio to first basis vector: {shortest_length/bkz_length:.3f}")
        
        except Exception as e:
            print(f"Enumeration failed or too slow: {e}")
    
    # Step 6: Summary
    print(f"\n6. SUMMARY")
    print("-" * 40)
    
    print(f"Lattice dimension: {dimension}")
    print(f"Total reduction time: {lll_time + bkz_time:.4f} seconds")
    print(f"")
    print(f"Quality progression:")
    print(f"  Initial:  ||b₁|| = {initial_length:.3f}, δ = {initial_hermite:.6f}")
    print(f"  LLL:      ||b₁|| = {lll_length:.3f}, δ = {lll_hermite:.6f}")
    print(f"  BKZ-{block_size}:   ||b₁|| = {bkz_length:.3f}, δ = {bkz_hermite:.6f}")
    print(f"")
    print(f"Theoretical comparison:")
    print(f"  LLL guarantee: δ ≤ (4/3)^(n-1)/4 = {(4/3)**((dimension-1)/4):.6f}")
    print(f"  Achieved: δ = {bkz_hermite:.6f}")
    print(f"  Gap to GH: {bkz_length/gh:.3f}")
    
    return A, M

# Run the complete analysis
print("Running complete lattice analysis workflow...")
print("This demonstrates all major fpylll capabilities in one example.")

try:
    lattice, gso = complete_lattice_analysis(dimension=20, bits=25)
    print(f"\n✓ Analysis completed successfully!")
    print(f"✓ Final lattice basis stored in 'lattice'")
    print(f"✓ GSO object stored in 'gso'")
    
except Exception as e:
    print(f"\n✗ Analysis failed: {e}")
    print("This might be due to dimension too large or other computational limits")

print(f"\n" + "="*60)
print("FPYLLL TUTORIAL COMPLETE")
print("="*60)
print("You now have a comprehensive understanding of:")
print("• Lattice theory and computational problems")
print("• fpylll library structure and capabilities") 
print("• LLL and BKZ reduction algorithms")
print("• Cryptanalytic applications")
print("• Performance optimization techniques")
print("")
print("Ready for advanced lattice-based research and applications!")