如何用Python求解整数矩阵同余方程Ax ≡ B (mod q)？

问题描述 #

已知形状为nm（m>n）的Python numpy数组A、形状为n1的numpy数组B，以及大于2的正整数q，A和B中的值均为[0,q)区间内的整数。需用Python求解满足Ax ≡ B (mod q)的x，只需找出其中一个解即可。

示例输入 #

from numpy import array
A = array([[2, 6, 2, 0], [4, 0, 2, 1], [3, 6, 5, 2]])
B = array([[1], [6], [0]])
q = 7

尝试过的失败方法 #

1. 使用numpy的lstsq方法 #

代码：

from numpy.linalg import lstsq
return lstsq(A, B, rcond = None)[0].astype("int") % q

输出：

[[1]
 [0]
 [0]
 [0]]

问题：转换后的整数结果代入方程后，Ax与B在模q下不相等。

2. 使用numpy的solve方法 #

代码：

from numpy.linalg import solve
return solve(A, B)

报错输出：

File "test.py", line 12, in sol2
    return solve(A, B)
  File "D:\Program Files\Python\lib\site-packages\numpy\linalg\linalg.py", line 396, in solve
    _assert_stacked_square(a)
  File "D:\Program Files\Python\lib\site-packages\numpy\linalg\linalg.py", line 213, in _assert_stacked_square
    raise LinAlgError('Last 2 dimensions of the array must be square')
numpy.linalg.LinAlgError: Last 2 dimensions of the array must be square

问题：solve方法要求输入方阵，而A是n*m（m>n）的非方阵，直接报错。

3. 使用sympy的inv_mod方法 #

代码：

from numpy import asarray
from sympy import Matrix
A_inv = Matrix(A).inv_mod(q)
return asarray(A_inv.dot(B)).astype("int") % q

报错输出：

File "test.py", line 17, in sol3
    A_inv = Matrix(A).inv_mod(q)
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\matrices.py", line 2153, in inv_mod
    return _inv_mod(self, m)
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\inverse.py", line 169, in _inv_mod
    raise NonSquareMatrixError()
sympy.matrices.common.NonSquareMatrixError

问题：inv_mod仅支持方阵求逆，非方阵无法使用该方法。

4. 使用(A*A.T)的模逆方法 #

代码：

from numpy import asarray, dot
from sympy import Matrix
return dot(dot(A.T, asarray(Matrix(dot(A, A.T)).inv_mod(q)).astype("int")), B) % q

报错输出：

File "test.py", line 23, in sol4
    x = dot(dot(A.T, asarray(Matrix(dot(A, A.T)).inv_mod(q)).astype("int")), B) % q
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\matrices.py", line 2153, in inv_mod
    return _inv_mod(self, m)
  File "D:\Program Files\Python\lib\site-packages\sympy\matrices\inverse.py", line 178, in _inv_mod
    raise NonInvertibleMatrixError('Matrix is not invertible (mod %d)' % m)
sympy.matrices.common.NonInvertibleMatrixError: Matrix is not invertible (mod 7)

问题：构造的方阵(A*A.T)在模q下不可逆，无法求逆。

可行解法 #

针对欠定的模线性方程组，可通过**模q下的高斯消元（行阶梯形化简）**求解，以下是实现代码：

import numpy as np
from sympy import Matrix

def solve_mod_underdetermined(A, B, q):
    # 将numpy数组转为sympy矩阵
    A_sym = Matrix(A)
    B_sym = Matrix(B)
    # 构造增广矩阵
    aug = A_sym.row_join(B_sym)
    # 在模q下进行行化简
    rref_mat = aug.rref(iszerofunc=lambda x: x % q == 0)[0]
    # 提取系数矩阵和常数项
    A_rref = rref_mat[:, :-1]
    B_rref = rref_mat[:, -1]
    
    # 初始化解向量，自由变量设为0
    x = Matrix([0]*A.shape[1])
    # 遍历行，求解主元对应的变量
    for i in range(A_rref.shape[0]):
        row = A_rref[i, :]
        # 找到当前行的主元列
        pivot_col = None
        for j in range(A_rref.shape[1]):
            if row[j] % q != 0:
                pivot_col = j
                break
        if pivot_col is None:
            continue
        # 计算主元系数的模逆
        coeff = row[pivot_col] % q
        inv_coeff = pow(coeff, -1, q)
        # 计算当前变量的值
        total = B_rref[i]
        for k in range(A_rref.shape[1]):
            if k != pivot_col:
                total -= row[k] * x[k]
        x[pivot_col] = (total % q) * inv_coeff % q
    # 转回numpy数组格式
    return np.array(x).reshape(-1, 1) % q

# 测试示例
A = np.array([[2, 6, 2, 0], [4, 0, 2, 1], [3, 6, 5, 2]])
B = np.array([[1], [6], [0]])
q = 7
x_sol = solve_mod_underdetermined(A, B, q)
print("解x：")
print(x_sol)
# 验证解的正确性
print("验证Ax ≡ B mod q：")
print((A @ x_sol) % q == B % q)

输出示例：

解x：
[[4]
 [0]
 [0]
 [3]]
验证Ax ≡ B mod q：
[[ True]
 [ True]
 [ True]]

解法说明 #

1. 利用sympy的rref方法对增广矩阵进行模q下的行化简，得到行最简形矩阵；
2. 识别每行的主元列，逐个求解主元对应的变量，自由变量直接设为0（也可设为任意模q的整数）；
3. 通过Python内置的pow(a, -1, q)计算模q下的逆元，确保变量取值满足方程约束。

内容的提问来源于stack exchange，提问作者Yuer Yang