sympy计算smith标准型

1.源码实现

import itertools
from sympy import *

def nchoosek(a, b, d=1, n=1):
    c = []
    for i in itertools.combinations(range(a,b+1,d),n):
        c.append(list(i))
    return c

def smith_form(C):
    #1---求lambda矩阵的行列式因子(开始)
    n = C.shape[0]
    h = [];

    for k in range(1, n+1):
        y = [];

        #从n个元素(第一个元素到第n个元素的所有元素)中取出k项，其中所有可能的组合放在p中(组合数)
        p = Matrix(nchoosek(1, n, n=k));
        q = p;
        m = p.shape[0]

        #print("p=")
        #print(p)

        #1-b将所有的i阶行列式求出
        for i in range(0, m):
            a = p.row(i)

            for j in range(0, m):
                b = q.row(j)
                B = []

                for s in range(0, k):
                    g = [];

                    for t in range(0, k):
                        u = a[s]
                        v = b[t]
                        g.append(C[u-1, v-1])

                    B.append(g)

                c = Matrix(B).det()

                if c != 0:
                    y.append(c);
        #1-b结束
        o = len(y)

        #1-a--求所有j阶行列式的最大公因子（开始）
        for j in range(0, o-1):
            y[j+1] = gcd(y[j],y[j+1]);
        #1-a--结束

        h.append(y[o-1]);
    #1---求lambda矩阵的行列式因子(结束)

    d = [];
    d.append(h[0])

    for i in range(1, n):
        d.append(expand(h[i]/h[i-1]));

    #2---将上面所求行列式因子进行因式分解，得到因式分解后的行列式因子
    k = len(d)

    for i in range(0, k):
        t = d[i]
        t = factor(t)
        d[i] = t

    return diag(d)


#声明变量x
x = symbols('x')

#输入lambda矩阵
C = Matrix([[-x+1, 2*x-1, x], [x, x**2, -x], [x**2+1, x**2+x-1, -x**2]])

#计算smith标准型
D = smith_form(C)

print(C)
print(D)

2.运行及结果

$ python3 example.py
Matrix([[-x + 1, 2*x - 1, x], [x, x**2, -x], [x**2 + 1, x**2 + x - 1, -x**2]])
Matrix([[1], [x], [-x*(x**2 + 1)]])