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)]])
网友评论