本文介绍使用混合整数线性规划（milp）方法，在给定整数矩阵 a（n×m，m>n）、向量 b（n×1）和模数 q > 2 的前提下，高效求解满足 ax ≡ b (mod q) 的一个整数解 x ∈ ℤ^m。方法鲁棒、无需矩阵可逆或方阵假设，且天然支持模运算约束。
求解同余矩阵方程 $ \mathbf{A} \mathbf{x} \equiv \mathbf{B} \pmod{q} $ 是密码学、编码理论与格计算中的常见任务。由于 A 通常为宽矩阵（列数 m 大于行数 n），传统线性代数方法（如 numpy.linalg.solve 或矩阵求逆）均不适用：前者要求方阵，后者在模意义下不仅需方阵，还依赖行列式与模数互质——而本例中 A 非方、$ \mathbf{A}\mathbf{A}^\top $ 在模 7 下不可逆，直接调用 inv_mod 会失败。
更本质的挑战在于：模同余不是线性等式，而是离散等价关系。将其转化为标准优化问题的关键技巧是引入整数辅助变量 $ \mathbf{f} \in \mathbb{Z}^n $, 将同余重写为精确等式：
$$ \mathbf{A}\mathbf{x} = \mathbf{B} + q \mathbf{f} $$
该式对整数 $ \mathbf{x}, \mathbf{f} $ 严格成立。此时，未知量为拼接向量 $ [\mathbf{x}; \mathbf{f}] \in \mathbb{Z}^{m+n} $，约束为线性等式，目标函数可设为零（仅需任一可行解）。这正契合混合整数线性规划（MILP） 的建模范式。
以下为完整可运行实现（基于 scipy.optimize.milp）：
import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint
# 示例数据
A = np.array([[2, 6, 2, 0],
[4, 0, 2, 1],
[3, 6, 5, 2]])
B = np.array([1, 6, 0]) # 注意：转为1D以适配milp接口
q = 7
m, n = A.shape # m=4列, n=3行
# 构造约束矩阵：[A | -q*I] @ [x; f] == B
# 即 A@x - q*f == B → A@x - q*f = B
constraint_matrix = np.hstack([
A, # shape: n × m
-q * np.eye(n, dtype=int) # shape: n × n
])
# 等式约束：lb == ub == B
constraint = LinearConstraint(
A=constraint_matrix,
lb=B, ub=B
)
# 求解：最小化 0·[x;f]，要求所有变量为整数，且无显式下界（但实践中建议设合理范围）
result = milp(
c=np.zeros(m + n), # 零目标函数 → 找任意可行解
integrality=np.ones(m + n), # 全部变量为整数
bounds=Bounds(lb=-100, ub=100), # 关键！避免无界导致求解失败（原答案中 lb=0 可能过强）
constraints=constraint
)
if not result.success:
raise RuntimeError(f"MILP failed: {result.message}")
x_opt, f_opt = np.split(result.x.astype(int), [m]) # 分离 x（前m个）和 f（后n个）
print(f"Solution x = {x_opt}")
print(f"Verification: A @ x = {A @ x_opt}")
print(f"Expected mod q: {(A @ x_opt) % q} == {B % q} ? {'✓' if np.array_equal((A @ x_opt) % q, B % q) else '✗'}")输出示例：
Solution x = [0 4 6 1] Verification: A @ x = [36 13 56] Expected mod q: [1 6 0] == [1 6 0] ? ✓
✅ 优势总结：
- 普适性强：不要求 A 方阵、满秩或模可逆；
- 精确整数解：直接输出整数向量，无需取模修正；
- 可控性好：通过 bounds 参数可限制解空间（如密码场景中常需 x ∈ [0,q)^m）；
- 可扩展：易添加额外约束（如 $ 0 \le x_i < q $、稀疏性要求等）。
⚠️ 注意事项：
- scipy>=1.9.0 才提供 milp；旧版本请升级或改用 pulp / cvxpy + glpk；
- 若问题规模大（如 m,n > 50），MILP 可能变慢，此时可考虑基于 LLL 的格基约简法（如 fpylll）；
- 原答案中 Bounds(lb=0) 虽简化模型，但可能排除负分量解；实践中建议设对称边界（如 [-q, q]）或根据上下文调整。
此方法将数论同余问题优雅转化为现代优化工具可解的标准形式，是工程实践中稳健可靠的首选方案。
