# Calculate the Smith normal form of an integer matrix

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 : 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!

• Sandbox entry was: codegolf.meta.stackexchange.com/questions/2140/… May 9 at 17:03
• May we take the dimensions of the matrix as input as well? May 9 at 17:20
• @chunes Yes, certainly. (I count that as part of the specification of the input matrix, especially for example in C where the usual representation would be a pointer to a flat array, and the dimensions to use to interpret it as a matrix.) May 9 at 17:22
• I was confused for a bit by your sample program, but the matrices returned by it are not the $P$ and $Q$ in the post. They instead satisfy $D = P A Q$. Separately, you might want to explain that $d_{11} | d_{22}$ means $d_{11}$ divides $d_{22}$ evenly as I don't think that notation is particularly common for people who don't study number theory, etc. May 9 at 21:20
• Oh, right, that's another possibility for the decomposition, that you list $P^{-1}, Q^{-1}$ instead of $P, Q$ - which would be the natural output for example if you keep an augmented matrix starting with identity matrix augmentations. It might also be the output that would be more interesting for some of the applications of Smith normal form decomposition. May 10 at 0:22

# 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\$$.

• Hmm, the divide by zero message at the link reminded me: there might be some issue when there are multiple zeros in the answer. For example, if I try to apply it to $\begin{bmatrix} 3&3&3 \\ 4&4&4 \\ 5&5&5 \end{bmatrix}$, it's showing an answer of {1, 0, Indeterminate}. May 10 at 15:19
• @DanielSchepler Not a problem with 2 more bytes... May 10 at 15:20
• Rest@#/Most@#/. 0/0->0&@Table[GCD@@Join@@#~Minors~k,{k,0,Tr[1^#]}]&. May 10 at 23:57
• 56
– att
May 14 at 5:43

# PARI/GP, 20 bytes

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

m->Vecrev(matsnf(m))


Attempt This Online!

• From what I can see, matsmf without the flag returns the list of diagonal entries though in the opposite order from my convention; and since I allow just the list of diagonal entries as an output format, all you would really need to do would be to reverse that list. May 10 at 0:17
• i.e. m->Vecrev(matsnf(m)). May 10 at 0:24

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

• It looks like I should be able to get a significant savings by using std::views::iota instead of the hand-rolled loops; unfortunately, I don't have a compiler version supporting that on my machine at the moment, and I haven't found online compilers with the Eigen library installed. May 10 at 19:16
• Some minor golfs: k=0 can be removed, and you can add the =0 to the initialization of k; the && in k<m&&k<n can be &; likewise (v<P||!P)&&v can be (v<P|!P)&v; and a set of {} can be removed at for(j=k+1;j<n;++j){if(M(i,j)%P) (near the end). There are probably a lot more saves to be made. Nice answer btw! May 12 at 7:31
• I'm not convinced on (v<P|!P)&v - if it happens P is zero then !P would evaluate to 1, so then if v happens to be nonzero but even... May 15 at 21:18
• v<P|!P&&v should be ok Nov 3 at 8:23

# R, 165163 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.

• -2 bytes by inlining the definition of g and renaming it to ?. May 14 at 19:31
• I think you can use det(t(m[i,j])) to save a couple bytes. Jun 1 at 14:53

# Octave, 237229227 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

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];
for i=1:(min(size(A)))
if d(i)~=0
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));


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