如何用Python实现500位精度的实矩阵求逆？

方案一：用mpmath库（最省心） #

mpmath是专门做任意精度数值计算的库，自带矩阵求逆功能，完全能满足500位精度要求，不需要自己造轮子。

• 操作步骤：
  1. 先把全局精度设为500位：mpmath.mp.dps = 500
  2. 把实矩阵转成mpmath的矩阵类型
  3. 直接调用mpmath.inverse()得到逆矩阵
  4. 验证的话，用原矩阵乘逆矩阵，再用mpmath.chop()去掉极小误差，看结果是否接近单位矩阵

代码示例：

import mpmath

# 设定500位精度
mpmath.mp.dps = 500

# 构造目标矩阵，可直接写数值或从文件读取
A = mpmath.matrix([
    [1.2345678901234567890123456789, 2.3456789012345678901234567890],
    [3.4567890123456789012345678901, 4.5678901234567890123456789012]
])

# 求逆
A_inv = mpmath.inverse(A)

# 验证结果：原矩阵乘逆矩阵应接近单位矩阵
print(mpmath.chop(A * A_inv))

方案二：用sympy做符号化高精度求逆 #

sympy的符号矩阵支持任意精度计算，适合需要精确表达式或高精度数值输出的场景：

• 操作步骤：
  1. 用sympy的Matrix构造矩阵，数值可使用有理数或高精度浮点数
  2. 调用矩阵的inv()方法求逆
  3. 用sympy.N(结果, 500)输出500位精度的数值结果

代码示例：

import sympy

# 用有理数构造矩阵，避免初始精度损失
A = sympy.Matrix([
    [sympy.Rational(12345678901234567890123456789, 10**27), sympy.Rational(23456789012345678901234567890, 10**27)],
    [sympy.Rational(34567890123456789012345678901, 10**27), sympy.Rational(45678901234567890123456789012, 10**27)]
])

# 求逆并输出500位精度
A_inv = A.inv()
print(sympy.N(A_inv, 500))

方案三：基于decimal手动实现高斯消元 #

如果不想依赖第三方库，可以用decimal结合高斯消元法自己写求逆逻辑，全程用Decimal类型计算，绝对不会出现转低精度float的问题：

• 核心思路：构造增广矩阵[A | I]（I为单位矩阵），通过行变换把左边的A转成单位矩阵，右边剩下的就是逆矩阵，所有运算都用Decimal类型。

代码示例（简化版，仅适用于可逆矩阵）：

from decimal import Decimal, getcontext

# 设定500位精度
getcontext().prec = 500

def matrix_inverse(matrix):
    n = len(matrix)
    # 构造增广矩阵：原矩阵拼接单位矩阵
    aug = [row + [Decimal(1) if i == j else Decimal(0) for j in range(n)] for i, row in enumerate(matrix)]
    
    for col in range(n):
        # 找主元行（当前列绝对值最大的行，避免数值不稳定）
        pivot_row = max(range(col, n), key=lambda r: abs(aug[r][col]))
        aug[col], aug[pivot_row] = aug[pivot_row], aug[col]
        
        # 主元所在行归一化
        pivot_val = aug[col][col]
        for j in range(col, 2*n):
            aug[col][j] /= pivot_val
        
        # 消去其他行的当前列元素
        for r in range(n):
            if r != col and aug[r][col] != 0:
                factor = aug[r][col]
                for j in range(col, 2*n):
                    aug[r][j] -= factor * aug[col][j]
    
    # 提取逆矩阵部分
    return [row[n:] for row in aug]

# 示例矩阵，全部用Decimal类型输入
A = [
    [Decimal('1.2345678901234567890123456789'), Decimal('2.3456789012345678901234567890')],
    [Decimal('3.4567890123456789012345678901'), Decimal('4.5678901234567890123456789012')]
]

A_inv = matrix_inverse(A)

# 验证：原矩阵乘逆矩阵应接近单位矩阵
def multiply_matrices(a, b):
    n = len(a)
    m = len(b[0])
    p = len(b)
    result = [[Decimal(0) for _ in range(m)] for _ in range(n)]
    for i in range(n):
        for k in range(p):
            if a[i][k] == 0:
                continue
            for j in range(m):
                result[i][j] += a[i][k] * b[k][j]
    return result

product = multiply_matrices(A, A_inv)
for row in product:
    # 打印时保留10位小数，应接近[1,0]和[0,1]
    print([round(x, 10) for x in row])

关于你提到的Stack Overflow帖子问题 #

那个帖子里的%运算符导致转float32的问题，是因为代码用了依赖普通float类型的操作，而上面三个方案全程使用高精度专属类型（mpmath.mpf、sympy符号数、Decimal），不会触发自动转低精度float的情况，完全规避了这个问题。

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