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.


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


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!

  • \$\begingroup\$ Sandbox entry was: codegolf.meta.stackexchange.com/questions/2140/… \$\endgroup\$ May 9, 2023 at 17:03
  • \$\begingroup\$ May we take the dimensions of the matrix as input as well? \$\endgroup\$
    – chunes
    May 9, 2023 at 17:20
  • \$\begingroup\$ @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.) \$\endgroup\$ May 9, 2023 at 17:22
  • 3
    \$\begingroup\$ 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. \$\endgroup\$ May 9, 2023 at 21:20
  • \$\begingroup\$ 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. \$\endgroup\$ May 10, 2023 at 0:22

6 Answers 6


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

-35 from att


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

  • \$\begingroup\$ 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}. \$\endgroup\$ May 10, 2023 at 15:19
  • \$\begingroup\$ @DanielSchepler Not a problem with 2 more bytes... \$\endgroup\$ May 10, 2023 at 15:20
  • 1
    \$\begingroup\$ Rest@#/Most@#/. 0/0->0&@Table[GCD@@Join@@#~Minors~k,{k,0,Tr[1^#]}]&. \$\endgroup\$
    – alephalpha
    May 10, 2023 at 23:57
  • \$\begingroup\$ 56 \$\endgroup\$
    – att
    May 14, 2023 at 5:43

PARI/GP, 20 bytes

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


Attempt This Online!


  • 2
    \$\begingroup\$ 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. \$\endgroup\$ May 10, 2023 at 0:17
  • 1
    \$\begingroup\$ i.e. m->Vecrev(matsnf(m)). \$\endgroup\$ May 10, 2023 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.)

#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;
    std::cout << "\nOutput:\n" << M << "\n\n";

int main() {
    Eigen::MatrixXi M(3,3);

    Eigen::MatrixXi M2(2,2);

    Eigen::MatrixXi M3(3,3);

    Eigen::MatrixXi M4(2,2);

    Eigen::MatrixXi M5(2,4);

    Eigen::MatrixXi M6(3,3);

    return 0;
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ May 10, 2023 at 19:16
  • \$\begingroup\$ 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! \$\endgroup\$ May 12, 2023 at 7:31
  • \$\begingroup\$ 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... \$\endgroup\$ May 15, 2023 at 21:18
  • \$\begingroup\$ v<P|!P&&v should be ok \$\endgroup\$
    – ceilingcat
    Nov 3, 2023 at 8:23

R, 165 163 159 bytes


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.


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.

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

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');
    % 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];

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

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

    [6,10; 10,15],
    [2 2  4 6 ;2 -2 0 -2],

SageMath, 26 24 bytes

lambda m:m.smith_form(0)

SageMath matrix as I/O.

Try it online.


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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.