求整数矩阵的Smith范式

Question

给定整数A的m \times n矩阵，存在m \times m矩阵P、m \times n矩阵D和n \times n矩阵Q，使得：

• A = P D Q。
• P和Q是单模矩阵(即可逆矩阵，其逆矩阵也是整数矩阵)；
• D是对角线；
• D的每个对角线项d_{ii}是非负的；
• d_{11} \mid d_{22} \mid \cdots \mid d_{nn} 。

此外，矩阵D在这种表示中是唯一的。

根据上下文的不同，A的Smith范式要么只是对角矩阵D，要么A的Smith范式分解是元组(P, D, Q)。对于这个挑战，你只需要计算D。

一种常用的计算D的方法是通过一种看起来像欧氏算法计算gcd和高斯消去的组合的算法--应用基本行和列运算，直到矩阵符合D的所需格式。另一种计算D的方法是，对于每个i，d_{11} d_{22} \cdots d_{ii}等于A的i\times i子矩阵(包括非连续子矩阵)的所有行列式的gcd。

挑战

您需要编写一个函数或程序来计算输入矩阵的Smith范式。输出可以是全矩阵D的形式，也可以是D的对角线条目列表的形式，在命令式编程语言中，也可以通过引用将矩阵写入函数并将矩阵修改为Smith范式。

规则

• 这是密码-高尔夫：最短代码获胜。
• 适用标准漏洞禁令。
• 如果这适用于您的语言，则不需要担心整数溢出。

示例

1 2 3       1 0 0
4 5 6  ->   0 3 0
7 8 9       0 0 0

6  10       1 0
10 15  ->   0 10

6 0  0        1 0  0
0 10 0   ->   0 30 0
0 0  15       0 0  30

2 2       2 0
2 -2  ->  0 4

2 2  4 6       2 0 0 0
2 -2 0 -2  ->  0 4 0 0

3 3 3      1 0 0
4 4 4  ->  0 0 0
5 5 5      0 0 0

注: Mathematica已经有一个内置的计算Smith范式的方法.因此，您可以使用它来处理测试用例：在网上试试！

Answer 1

Wolfram语言(数学)，没有SNF内置，56个字节

-35来自att

(i=0;(!i+1)/1~Max~!i++&/.!a_:>GCD@@Join@@#~Minors~a)/@#&

在网上试试！

来自维基百科：

不变因子\alpha_i可以计算为\alpha_i = \frac{d_i(A)}{d_{i-1}(A)}，其中d_i(A)等于矩阵A和d_0(A)=1的所有i\times i次元的行列式的最大公因子。

Answer 2

帕里/GP，20字节

-25字节，这要归功于@丹尼尔·谢普勒和@Parcly Taxel。

m->Vecrev(matsnf(m))

在网上试试！

相关。

Answer 3

C++和特征3库，470个字节

由于我们没有收到使用高斯消除风格的解决方案提交，我想我会在实现和高尔夫球解决方案自己。(我不会把一个自我回答算作竞争对手，但如果有人想把它捡起来，然后用它跑，那我就没问题了。)

#include<Eigen/Core>
#define r M.row
#define c M.col
#define G goto R;
void f(auto&M){int m=M.rows(),n=M.cols(),i,j,k,P,v;for(k=0;k<m&&k<n;++k){R:if(M(k,k)<0)r(k)*=-1;P=M(k,k);for(i=k;i<m;++i)for(j=k;j<n;++j){v=abs(M(i,j));if((v<P||!P)&&v){r(k).swap(r(i));c(k).swap(c(j));G}}if(P){for(i=k+1;i<m;++i){r(i)-=M(i,k)/P*r(k);if(M(i,k))G}for(j=k+1;j<n;++j){c(j)-=M(k,j)/P*c(k);if(M(k,j))G}for(i=k+1;i<m;++i)for(j=k+1;j<n;++j){if(M(i,j)%P){r(i)+=r(k);c(j)-=M(i,j)/P*c(k);G}}}}}

对于某些测试代码，请追加以下内容：

#include <iostream>
void test(Eigen::MatrixXi M) {
    std::cout << "Input:\n" << M;
    f(M);
    std::cout << "\nOutput:\n" << M << "\n\n";
}

int main() {
    Eigen::MatrixXi M(3,3);
    M<<1,2,3,4,5,6,7,8,9;
    test(M);

    Eigen::MatrixXi M2(2,2);
    M2<<6,10,10,15;
    test(M2);

    Eigen::MatrixXi M3(3,3);
    M3<<6,0,0,0,10,0,0,0,15;
    test(M3);

    Eigen::MatrixXi M4(2,2);
    M4<<2,2,2,-2;
    test(M4);

    Eigen::MatrixXi M5(2,4);
    M5<<2,2,4,6,2,-2,0,-2;
    test(M5);

    Eigen::MatrixXi M6(3,3);
    M6<<3,3,3,4,4,4,5,5,5;
    test(M6);

    return 0;
}