Calculate the Smith normal form of an integer matrix

Question

Given an \$m \times n\$ matrix of integers A, there exist a \$m \times m\$ matrix P, an \$m \times n\$ matrix D, and an \$n \times n\$ matrix Q such that:

• \$A = P D Q\$.
• P and Q are unimodular matrices (i.e. matrices which are invertible and whose inverses are also integer matrices);
• D is diagonal;
• each diagonal entry \$d_{ii}\$ of D is nonnegative; and
• \$d_{11} \mid d_{22} \mid \cdots \mid d_{nn} \$.

Furthermore, the matrix D is unique in this representation.

Depending on context, the Smith normal form of A is either just the diagonal matrix D, or a Smith normal form decomposition of A is the tuple \$(P, D, Q)\$. For this challenge, you only need to calculate D.

One common way to calculate D is via an algorithm that looks like a combination of the Euclidean algorithm for calculating gcd and Gaussian elimination -- applying elementary row and column operations until the matrix is in the desired format for D. Another way to calculate D is to use the fact that for each i, \$d_{11} d_{22} \cdots d_{ii}\$ is equal to the gcd of all determinants of \$i\times i\$ submatrices (including non-contiguous submatrices) of A.

The challenge

You are to write a function or program that calculates the Smith normal form of an input matrix. The output may either be in the form of the full matrix D, or in the form of a list of the diagonal entries of D. In an imperative programming language, it is also acceptable to write a function that takes a matrix by reference and mutates the matrix into its Smith normal form.

Rules

• This is code-golf: shortest code wins.
• Standard loophole prohibitions apply.
• You do not need to worry about integer overflow, if that applies in your language.

Examples

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

Note: Mathematica already has a built-in to calculate Smith normal form. As such, you can use this to play around with test cases: Try it online!

Answer 1

Wolfram Language (Mathematica), no SNF builtin, 56 bytes

-35 from att

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

Try it online!

From Wikipedia:

The invariant factors \$\alpha_i\$ can be computed as $$\alpha_i = \frac{d_i(A)}{d_{i-1}(A)}$$ where \$d_i(A)\$ equals the greatest common divisor of the determinants of all \$i\times i\$ minors of the matrix A and \$d_0(A)=1\$.

Answer 2

PARI/GP, 20 bytes

-25 bytes thanks to @Daniel Schepler and @Parcly Taxel.

m->Vecrev(matsnf(m))

Attempt This Online!

Related.

Answer 3

C++ with the Eigen3 library, 470 bytes

Since we were not getting submissions using the Gaussian elimination style solution, I thought I would take a crack at implementing and golfing a solution myself. (I won't count a self-answer as a competitor, but if somebody wants to pick this up and run with it, that would be fine with me.)

#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}}}}}

For some test code, append this:

#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;
}

Answer 4

R, 165 163 159 bytes

\(m,x,y,R=Reduce,C=combn)R(\(l,k)R(`?`<-\(x,y,z=y%%x)`if`(x&z,z?x,x)+y*!x,sort(abs(C(x,k,\(i)C(y,k,\(j)round(det(t(m[i,j]))))))))/max(1,l),1:min(x,y),1,,T)[-1]

Attempt This Online!

-2 bytes thanks to pajonk and further -4 thanks to Giuseppe.

Takes input as matrix m and its dimensions x,y, outputs SNF diagonal.

Explanation

Uses the same algorithm as Parcly Taxel's answer.

For each k in 1:min(x,y) compute the minors:

combn(x,k,\(i)...) generate all combinations i of k row indices

combn(y,k,\(j)...) generate all combinations j of k column indices

round(det(m[i,j])) calculate the determinants of submatrices and round, as otherwise floating point inaccuracies lead to problems in GCD

sort(abs(...) discard the signs and sort in order expected by GCD

Reduce('?'<-..., ...) Reduce by GCD. Overall, a lot of bytes are wasted due to no built-in GCD in R...

Reduce(\(l,k)) .../max(1,l), init=1, accumulate=T) Divide each term by the previous one, unless it's zero

..[-1] Drop the initializer 1 and return.

Answer 5

Octave, 237 229 227 225 bytes.

Saved 10 bytes thanks to the comment of @ceilingcat

---

Use the formula in Wikipedia.

The invariant factors \$\alpha_i\$ can be computed as $$\alpha_i = \frac{d_i(A)}{d_{i-1}(A)} $$ where \$d_i(A)\$ equals the greatest common divisor of the determinants of all \$i\times i\$ minors of the matrix A and \$d_0(A)=1\$.
\$ \text{Note that here } d_i(A)\geq 1 \$

---

Golfed version. Try it online!

function a=z(A)[m,n]=size(a=A*0);d=1;for i=1:min(m,n);q=[];for r=combnk(1:m,i)';for c=combnk(1:n,i)';q=[q;det(A(r,c))];end;end;for j=1:numel(q);x=abs(round(q(j)));r={x gcd(r,x)}{(j>1)+1};end;if d~=0;a(i,i)=r/d;end;d=r;end;end

Ungolfed version. Try it online!

clear all;close all;clc;

function dets = minor_det(A, i)
    % Check if A is a square matrix or not
    [m, n] = size(A);
    
    % If i is greater than any dimension of A, return an error
    if i > m || i > n
        error('i cannot be greater than the dimensions of A');
    end
    
    % Get all combinations of i rows and columns
    rows = combnk(1:m, i);
    cols = combnk(1:n, i);
    
    % Initialize the array of determinants
    dets = [];
    
    % For each combination of rows and columns, compute the determinant
    for r = 1:size(rows, 1)
        for c = 1:size(cols, 1)
            minor = A(rows(r,:), cols(c,:));
            det_minor = det(minor);
            dets = [dets; det_minor];
        end
    end
end

function answer = func(A)
d=[];
for  i=1:min(size(A))
    minor_determinants = minor_det(A,i);
    res = round(abs(minor_determinants(1)));
    for index = 2:numel(minor_determinants)
        res = gcd(res, round(abs(minor_determinants(index))));
    end
    d(end+1)=res;
end
d=[1,d];
answer=zeros(size(A));
for i=1:(min(size(A)))
    if d(i)~=0
        answer(i,i)=d(i+1)/d(i);
    end
end
end

function out = map_func(matrices)
out = cellfun(@(matrix) func(matrix), matrices, 'UniformOutput', false);
end

matrices={
    [1,2,3;4,5,6;7,8,9],
    [6,10; 10,15],
    [6,0,0;0,10,0;0,0,15],
    [2,2;2,-2],
    [2 2  4 6 ;2 -2 0 -2],
    [3,3,3;4,4,4;5,5,5]
    };
disp(map_func(matrices));

Answer 6

SageMath, 26 24 bytes

lambda m:m.smith_form(0)

SageMath matrix as I/O.

Try it online.

Explanation:

Should speak for itself, but I'll add it anyway.

lambda m:          # Lambda function with matrix as argument:
  m.smith_form( )  #  Get the Smith Normal Form of the given matrix,
               0   #  but only return its diagonal matrix D

smith_form() will output a tuple of three matrices \$(D,P,Q)\$ by default. But by overriding its first argument with False, it'll result in just its diagonal matrix \$D\$, hence the need for the 0-argument.