从零实现LLL算法：用Python破解格密码与背包问题

在密码学与数学的交汇处，格(lattice)理论正掀起一场静默的革命。当量子计算机威胁着传统加密体系时，基于格的密码方案因其抗量子特性而备受瞩目。但今天我们不谈抽象理论，而是打开Python编辑器，亲手实现格基约化的核心武器——LLL算法，让它帮我们完成三项实战任务：分解多项式、发现整数关系，以及破解简易背包密码系统。

1. 环境配置与数学工具准备

工欲善其事，必先利其器。我们先搭建实验环境，安装必要的数学库：

pip install numpy sympy matplotlib

这三个库将构成我们的计算核心：

• NumPy ：处理高精度数值计算和矩阵运算
• SymPy ：进行符号计算和多项式操作
• Matplotlib ：可视化算法中间过程

理解格的基本概念至关重要。一个n维格可以看作是一组线性无关向量（称为格基）的所有整数线性组合构成的集合。用数学表达式表示就是：

$$ L = { \sum_{i=1}^{n} a_i \mathbf{b}_i \mid a_i \in \mathbb{Z} } $$

其中$\mathbf{b}_1, ..., \mathbf{b}_n$是格基向量。LLL算法的核心目标就是找到一组"优质"的基向量——它们相对短且接近正交。

2. LLL算法核心实现

让我们从Gram-Schmidt正交化开始，这是LLL的基础步骤。以下函数计算一组向量的正交基：

def gram_schmidt(basis):
    ortho_basis = []
    for i in range(len(basis)):
        v = basis[i].copy()
        for j in range(i):
            mu = np.dot(basis[i], ortho_basis[j]) / np.dot(ortho_basis[j], ortho_basis[j])
            v -= mu * ortho_basis[j]
        ortho_basis.append(v)
    return ortho_basis

完整的LLL算法实现需要处理两个关键条件：大小缩减条件和Lovász条件。以下是Python实现的核心部分：

def LLL(B, delta=0.75):
    n = len(B)
    B = np.array(B, dtype=float)
    GS = gram_schmidt(B)
    
    k = 1
    while k < n:
        # 大小缩减步骤
        for j in range(k-1, -1, -1):
            mu = np.dot(B[k], GS[j]) / np.dot(GS[j], GS[j])
            if abs(mu) > 0.5:
                B[k] -= round(mu) * B[j]
                GS = gram_schmidt(B)
        
        # Lovász条件检查
        norm_GSk = np.dot(GS[k], GS[k])
        norm_GSk1 = np.dot(GS[k-1], GS[k-1])
        mu_k_k1 = np.dot(B[k], GS[k-1]) / np.dot(GS[k-1], GS[k-1])
        
        if norm_GSk >= (delta - mu_k_k1**2) * norm_GSk1:
            k += 1
        else:
            B[[k-1, k]] = B[[k, k-1]]
            GS = gram_schmidt(B)
            k = max(k-1, 1)
    
    return B

参数 delta (δ)控制约化质量，通常取值在(0.25,1)之间。δ越接近1，得到的基质量越高，但计算成本也越大。

3. 实战应用一：多项式分解

考虑分解多项式$f(x) = x^6 + 4x^5 - 12x^4 + 10x^3 - 7x^2 + 2x - 1$。我们可以构造特定格，用LLL算法寻找其因式：

from sympy import symbols, Poly

def factor_poly_with_LLL(f):
    x = symbols('x')
    degree = f.degree()
    # 构造格矩阵
    lattice = []
    for i in range(degree):
        row = [0]*i + [1] + [0]*(degree-1-i)
        row.append(f.coeff_monomial(x**i))
        lattice.append(row)
    lattice.append([0]*degree + [1])
    
    # 应用LLL
    reduced_basis = LLL(np.array(lattice))
    
    # 寻找可能的因式
    for vec in reduced_basis:
        if vec[-1] == 0:
            g = sum(c*x**i for i,c in enumerate(vec[:-1]))
            if g != 0 and f % g == 0:
                return g, f/g
    return f, 1

运行这个函数，我们可能发现$f(x) = (x^3 + 2x^2 - 3x + 1)(x^3 + 2x^2 - x - 1)$。这种技术在密码分析中非常有用，特别是当处理某些特殊形式的密钥时。

4. 实战应用二：发现整数关系

Machin公式$\pi/4 = 4\arctan(1/5) - \arctan(1/239)$展示了优美的整数关系。让我们用LLL重新发现它：

def find_integer_relation(numbers, precision=50):
    from decimal import Decimal, getcontext
    getcontext().prec = precision
    
    scaled = [int(Decimal(str(n)) * 10**precision) for n in numbers]
    n = len(scaled)
    
    # 构造格矩阵
    lattice = np.eye(n+1, dtype=int)
    for i in range(n):
        lattice[i,-1] = scaled[i]
    
    # 应用LLL
    reduced = LLL(lattice)
    
    # 寻找最小向量
    min_vec = min(reduced, key=lambda v: sum(x*x for x in v[:-1]))
    return min_vec[:-1]

输入$\arctan(1/5)$和$\arctan(1/239)$的系数，算法可能返回[1, -4, 1]，对应关系$4a - b - c = 0$，这正是Machin公式的变体。

5. 实战应用三：破解背包密码

背包密码系统基于子集和问题的难解性。给定公钥(超递增序列的线性组合)$\mathbf{a}$和密文$S = \sum a_i x_i$，其中$x_i \in {0,1}$，攻击步骤如下：

1. 构造格基矩阵： $$ B = \begin{bmatrix} 1 & 0 & \cdots & 0 & Na_1 \ 0 & 1 & \cdots & 0 & Na_2 \ \vdots & \vdots & \ddots & \vdots & \vdots \ 0 & 0 & \cdots & 1 & Na_n \ 0 & 0 & \cdots & 0 & -NS \end{bmatrix} $$

2. 应用LLL算法寻找短向量

3. 从短向量中恢复二进制解

Python实现：

def attack_knapsack(public_key, S, N=100):
    n = len(public_key)
    basis = np.eye(n+1)
    for i in range(n):
        basis[i, -1] = N * public_key[i]
    basis[-1, -1] = -N * S
    
    reduced = LLL(basis)
    
    for vec in reduced:
        if set(vec[:-1]).issubset({0,1,-1}) and vec[-1] == 0:
            return [abs(x) for x in vec[:-1]]
    return None

例如，对于公钥[183, 241, 87, 94, 250]和密文548，算法可能返回解[1,1,0,1,0]，对应明文11010。

6. 算法优化与可视化分析

原始LLL算法有多个优化方向。我们可以添加进度可视化来观察算法运行过程：

def visualize_LLL(B, delta=0.75):
    history = []
    
    def callback(k, B_curr):
        history.append((k, B_curr.copy()))
    
    # 修改后的LLL实现会调用callback
    final_basis = LLL_with_callback(B, delta, callback)
    
    # 绘制向量长度变化
    plt.figure(figsize=(10,6))
    for i in range(B.shape[1]):
        lengths = [np.linalg.norm(step[1][:,i]) for step in history]
        plt.plot(lengths, label=f'Vector {i+1}')
    plt.xlabel('Iteration')
    plt.ylabel('Vector Length')
    plt.legend()
    plt.show()
    return final_basis

观察向量长度变化可以帮助理解算法如何逐步改进基的质量。典型情况下，我们会看到向量长度震荡下降，最终稳定在较短的长度。

另一个重要优化是使用浮点运算的数值稳定性处理。原始实现可能遇到精度问题，特别是处理高维格时。我们可以引入以下改进：

def improved_LLL(B, delta=0.75, precision=1e-10):
    B = np.array(B, dtype=np.float64)
    n = B.shape[0]
    GS = gram_schmidt(B)
    
    k = 1
    while k < n:
        # 改进的大小缩减，处理数值误差
        for j in range(k-1, -1, -1):
            mu = np.dot(B[k], GS[j]) / (np.dot(GS[j], GS[j]) + precision)
            if abs(mu) > 0.5 + precision:
                B[k] -= np.round(mu) * B[j]
                GS = gram_schmidt(B)
        
        # 稳定的Lovász条件检查
        norm_k = np.dot(GS[k], GS[k]) + precision
        norm_k1 = np.dot(GS[k-1], GS[k-1]) + precision
        mu = np.dot(B[k], GS[k-1]) / norm_k1
        
        if norm_k >= (delta - mu**2) * norm_k1 - precision:
            k += 1
        else:
            B[[k-1, k]] = B[[k, k-1]]
            GS = gram_schmidt(B)
            k = max(k-1, 1)
    
    return B

这种改进虽然增加了少量计算开销，但显著提高了算法在高维情况下的可靠性。