# (RGS 4/5) Inverting matrices modulo m

Given an integer matrix M and a modulus m, find an inverse of M modulo m. If the matrix M is not invertible modulo m, the behaviour is left unspecified.

# Matrix inverse

If M is a square matrix, its inverse exists if and only if its determinant is not 0. Similarly, when we are talking about matrices modulo m, the inverse of M will exist if and only of the determinant of M is invertible modulo m, which happens when the determinant is coprime with m.

The inverse of M is a square matrix inv(M) such that M*inv(M) = inv(M)*M = I, where

$$I = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 1 \end{bmatrix}$$

has the same shape of M and is called the identity matrix. As an example, consider the first test case, where [[22, 43], [29, 37]] is the inverse of [[26, 16], [38, 41]] mod 45:

$$\begin{bmatrix} 26&16\\38&41 \end{bmatrix} \cdot \begin{bmatrix} 22&43\\29&37 \end{bmatrix} = \begin{bmatrix} 1036&1710\\2025&3151 \end{bmatrix} \equiv \begin{bmatrix} 1 & 0 \\ 0&1 \end{bmatrix} \mod 45$$

# Input

A square matrix M with integer values between 0 and m-1, inclusive, and a positive integer m > 1. The matrix may be given in any sensible format, including

• a list of lists, where the inner lists encode the rows, like M = [[1, 2], [3, 4]], or a flattened version, like M = [1, 2, 3, 4]
• a list of lists, where the inner lists encode the columns, like M = [[1, 3], [2, 4]], or a flattened version, like M = [1, 3, 2, 4]

where these encode the matrix

$$\\$$$$\begin{bmatrix} 1 & 2 \\ 3 & 4 \end{bmatrix}$$$$\\$$

The integer m giving the modulus.

You may also accept the size of the matrix as input.

The inputs can be given in any order.

# Output

A matrix representing the inverse of M modulo m. You may assume such an inverse exists. Preferable format is for each matrix entry $$\a_{i,j}\$$ to satisfy $$\0 \leq a_{i,j} < m\$$ but this is just to make it easier to compare with the test cases.

# Test cases

45, [[26, 16], [38, 41]] -> [[22, 43], [29, 37]]

39, [[29, 50], [29, 1]] -> [[16, 19], [4, 35]]

35, [[24, 14], [48, 45]] -> [[5, 7], [4, 33]]

53, [[43, 20], [15, 8]] -> [[5, 14], [37, 7]]

49, [[15, 11, 30], [20, 12, 40], [33, 25, 2]] -> [[33, 28, 23], [25, 18, 0], [25, 48, 13]]

37, [[8, 9, 22, 17], [24, 30, 30, 19], [39, 8, 45, 23], [5, 30, 22, 33]] -> [[18, 17, 26, 20], [29, 36, 23, 1], [19, 0, 9, 3], [30, 23, 14, 21]]


This is so shortest submission in bytes, wins! If you liked this challenge, consider upvoting it... And happy golfing!

This is the fourth challenge of the RGS Golfing Showdown. If you want to participate in the competition, you have 96 hours to submit your eligible answers. Remember there is still 300 reputation in prizes! (See 6 of the rules)

Also, as per section 4 of the rules in the linked meta post, the "restricted languages" for this third challenge are only Jelly, V (vim) and 05AB1E so submissions in these languages are not eligible for the final prize. But they can still be posted!!

Otherwise, this is still a regular challenge, so enjoy!

• I don't think it's quite true that a matrix is invertible if it has a nonzero determinant, when working modulo m where m is composite. For example, [[2, 0], [0, 1]] has determinant 2 but isn't invertible modulo 6.
– xnor
Mar 7 '20 at 8:21
• @xnor fixed, what I really wanted was "coprime with m".
– RGS
Mar 7 '20 at 8:23
• May we take the size of the matrix as input? Mar 7 '20 at 9:00
• In the description, you say A square matrix M with integer values between 0 and m-1, inclusive. However, in the last test case, one of the elements of the matrix is 39, greater than m, which is 37 Mar 8 '20 at 5:35
• It might be nice to have a test case where at least the top left element is not relatively prime to the modulus - to force failure in a naive Gaussian elimination. In fact, it would be nice to have a test case where all elements of the original matrix are not relatively prime to the modulus - for example, modulus 6, matrix [[2, 3], [3, 2]] Mar 9 '20 at 19:03

# R, 68 bytes

function(M,m,n,A=M){while(any(A%*%M%%m!=diag(n)))A[]=rpois(n^2,9)
A}


Try it online!

Strikingly slow. Will most likely time out for all test cases on TIO, but is guaranteed to give an answer eventually.

Works by rejection sampling: generates random matrices A, with each value taken from a $$\Poisson(9)\$$ distribution, until a solution is found.

Note that to get A of the correct dimensions, it is 6 bytes shorter to initialize it as A=M and then replace all the values with A[]=rpois(n^2,9) than to create it directly with A=matrix(rpois(n^2,9),n).

• small typo, this should return A rather than I, right? Assuming it ever returns, haha. Mar 10 '20 at 19:09
• @Giuseppe Of course, it should return A. Thanks! Mar 10 '20 at 19:26

# J, 18 16 bytes

(]%1+.]^5 p:[)%.


Try it online!

Resolves p/q mod n element-wise (instead of using det(M) to resolve the modular inverse globally).

Abuses GCD of rational numbers to extract 1/q from p/q.

### How it works

(]%1+.]^5 p:[)%.    NB. left arg = modulo, right arg = matrix
(            )%.    NB. bind inv(matrix) as new right arg
5 p:[       NB. phi(modulo)
]^            NB. inv(matrix)^phi(modulo) element-wise
1+.              NB. GCD with 1; GCD(1, p/q) = 1/q
]%                 NB. Divide inv(matrix) by the above element-wise


# J, 18 bytes

%.@]*-/ .*@]^5 p:[


Try it online!

A dyadic tacit function that takes modulo (left arg) and the matrix (right arg), and gives possibly very large-valued modular inverse of the matrix. To reduce the range, prepend [| at the start of the function.

### How it works: the math

A simple mathematical way to calculate the modular inverse of a matrix is the following:

\begin{align} M^{-1} \text{ mod }n &= \text{cofactor}(M) \times \bigl((\det M)^{-1} \text{ mod }n \bigr) \\ &= M^{-1} \times \det M \times \bigl((\det M)^{-1} \text{ mod }n \bigr) \end{align}

If the matrix $$\M\$$ is invertible modulo $$\n\$$, we know that $$\(\det M)^{-1} \text{ mod }n\$$ exists, and it can be found using Euler's theorem:

$$(\det M)^{-1} \equiv (\det M)^{\varphi(n)-1} \text{ mod }n$$

Then we can simplify the original equation to

\begin{align} M^{-1} \text{ mod }n &= M^{-1} \times \det M \times \bigl((\det M)^{\varphi(n)-1} \text{ mod }n \bigr) \\ &\equiv M^{-1} \times (\det M)^{\varphi(n)} \mod{n} \end{align}

And now the fun fact: J has built-ins for matrix inverse, matrix determinant, and Euler's totient function. And it uses built-in rational numbers when calculating the matrix inverse!

### How it works: the code

%.@]*-/ .*@]^5 p:[    NB. left arg = modulo, right arg = matrix
5 p:[    NB. totient(modulo)
-/ .*@]          NB. det(matrix)
^         NB. det(matrix) ^ totient(modulo)
%.@]                  NB. inv(matrix)
*                 NB. inv(matrix) * det(matrix) ^ totient(modulo)

• Nice, I knew there had to be some way to get modular matrix inversion from regular matrix inversion.
– xnor
Mar 10 '20 at 11:25
• I'm wondering if there's a way to do this without computing the matrix determinant. If I understand what's going on, the matrix inverse gives rational entries like p/q, and we need to re-interpret each of them as p/q done modulo n producing integer values. It would suffice to multiply everything by any large value that clears all the denominators and thereby produces integers, and equals 1 modulo n. One such value is indeed det(M)^phi(n). But maybe there's a shorter one.
– xnor
Mar 10 '20 at 11:46
• @xnor Your understanding is correct. Without det(M), the only other way would be to compute p/q mod n element-wise (or entirely fall back to Gaussian elimination). Mar 10 '20 at 13:58

# Wolfram Language (Mathematica), 23 bytes

¯\_(ツ)_/¯ the answer was in the documentation of Modulus

Inverse[#2,Modulus->#]&


Try it online!

# JavaScript (ES6),  209  206 bytes

Takes input as (modulo)(matrix).

This transposes the matrix of cofactors (resulting in the adjugate) and multiply it by the inverse of the determinant of $$\M\$$ modulo $$\m\$$.

m=>M=>M.map((r,y)=>r.map((_,x)=>((g=k=>(++k*D(M)%m+m)%m-1?g(k):x+y&1?-k:k)*D(h(M,x).map(r=>h(r,y)))%m+m)%m),h=(a,n)=>a.filter(_=>n--),D=M=>+M||M.reduce((s,[v],i)=>s+(i&1?-v:v)*D(h(M,i).map(r=>h(r,0))),0))


Try it online!

## Commented

### Helper function $$\h\$$

The function $$\h\$$ removes the $$\n\$$-th entry from the array $$\a[\:]\$$.

h = (a, n) =>                // a[] = array, n = index
a.filter(_ => n--)         // keep all but the n-th entry


### Helper function $$\D\$$

The function $$\D\$$ computes the determinant of the matrix $$\M\$$.

D = M =>                     // M[] = input matrix
+M ||                      // if M[] is 1x1, stop recursion and return its unique value
M.reduce((s, [v], i) =>    // otherwise, for each value v at (0, i):
s +                      //   add to the sum
(i & 1 ? - v : v) *      //   either v or -v depending on the parity of i
D(                       //   multiplied by the result of a recursive call with:
h(M, i)                //     M[] without the i-th row
.map(r => h(r, 0))     //     and without the first column
),                       //   end of recursive call
)                          // end of reduce()


### Main function

m => M =>                    // m = modulo, M[] = matrix
M.map((r, y) =>            // for each position y:
r.map((_, x) =>          //   for each position x:
(                      //
( g = k =>           //     g is a recursive function taking a counter k
( ++k *          //       increment k and multiply it
D(M)           //       by the determinant of M
% m + m        //
) % m - 1 ?      //       if it's not congruent to 1 modulo m:
g(k)           //         try again until it is
:                //       else:
x + y & 1 ? -k //         return either k or -k
: k  //         depending on the parity of x+y
) *                //     initial call to g with a zero'ish value
D(                   //     multiply by the determinant of:
h(M, x)            //       M[] without the x-th row
.map(r => h(r, y)) //       and without the y-th column
) % m + m            //     return the result modulo m
) % m                  //
)                        //   end of inner map()
)                          // end of outer map()


# Jelly, 25 bytes

ÆḊ×Ɱ⁹%ỊTḢ×ZÆḊ-Ƥ$-ƤNÐe⁺€Zʋ  Try it online! A dyadic link taking the matrix as its left argument and the modulus as its right. Returns a matrix. Append a % to get it within the range 0, m # SageMath, 48 33 bytes Saved 15 bytes thanks to ovs!!! lambda m,M:~Matrix(Integers(m),M)  Nothing on TIO for SageMath unfortunately. Modular inverse of a matrix M (input as a Python list of lists) mod m. • @ovs So sorry, is that correct? On my phone so can't check. Mar 7 '20 at 14:43 • Yes this works. This is the same operator as you use in expressions like -~n. – ovs Mar 7 '20 at 14:46 • @ovs Makes sense - thanks! :-) Mar 7 '20 at 14:52 • Is this a library for Python or another language? Mar 7 '20 at 16:48 • @S.S.Anne If you follow the above link you'll get to the dev page. It's basically a whole system built on top of Python w numpy/scipy that does symbolic math!! Mar 7 '20 at 17:06 # Sledgehammer, 6 bytes ⠑⡿⡆⠱⣁⣭  Decompresses into this Wolfram Language function: Inverse[#2, Modulus -> #1]  Try it online! # Charcoal, 41 bytes ＦＥＸθ×ηη⪪Ｅ×ηη÷ιＸθλη¿⬤ι⬤ζ⁼⁼λν﹪ΣＥμ×ξ§§ιπλθＩι  Try it online! Link is to verbose version of code. Takes input as $$\ m, n, M \$$ where $$\ n \$$ is the size of $$\ M \$$, and does not reduce its output modulo $$\ m \$$ (can be done at a cost of 2 bytes). Stupidly slow, so don't try this with realistic values. Explanation: ＦＥＸθ×ηη⪪Ｅ×ηη÷ιＸθλη  There are $$\ m^{n^2} \$$ possible square matrices of size $$\ n \$$ with coefficients between $$\ 0 \$$ and $$\ m \$$. Looping over this value, compute each matrix, but don't bother reducing the terms modulo $$\ m \$$. Then, loop over the list of matrixes. ¿⬤ι⬤ζ⁼⁼λν﹪ΣＥμ×ξ§§ιπλθ  Perform the steps of matrix multiplication of this matrix by the input matrix, reduce it modulo $$\ m \$$, and compare the each result to the appropriate value of the identity matrix. Ｉι  If this was the inverse then print the matrix. • By stupidly slow, I mean it can only cope with 4x4 with a modulus of 2 or 3x3 with a modulus of 2 or 3 or 2x2 with a modulus of up to 21 (23 at a cost of 2 bytes) on TIO. – Neil Mar 7 '20 at 22:42 # MATL, (25?) 31 29 26 bytes My first MATL answer -5 bytes & a bug-fix (+2) thanks to Luis Mendo! The trailing . may be unnecessary - it is if there is only ever a single inverse of M with elements modulo m. :inZ^!"&G@[]eY*w\tZyXy=?@.  A full program which prints the elements in row major order separated by newlines. Try it online! - Too slow for any of the given test cases. Quite possibly not the best approach for MATL. ### How? :inZ^!"&G@[]eY*w\tZyXy=?@. - expects inputs m and M : - range (m) -> [1,2,...,m] i - push input (M) n - number of elements Z^ - ([1,2,...,m]) Cartesian power (#elements(M)) ! - transpose " - for each column, C: &G - push both inputs @ - push C [] - push an empty array (to make e work as below) e - reshape (C) to square matrix of side ceil(#elements(C)^0.5) Y* - (reshaped C) matrix multiplication (copy of M) w - swap top two stack entries \ - (multiplication result) modulo (copy of m) t - duplicate top of stack Zy - size Xy - (size by size) identity matrix = - equal -> logical matrix ? - if all are truthy: @ - push C . - break - implicit print of stack (the valid C)  • @LuisMendo - I'm almost certain you'll beat this using another method, but any tips for this one would be welcome! Mar 7 '20 at 17:48 • ...also why is w "necessary" where I have it / is there something more normal that should go there? Mar 7 '20 at 17:52 • Thanks! Yes the 2 was somehow lost and should have been the column size. You're also correct about q being unnecessary since we don't have to give an answer modulo m. Mar 10 '20 at 20:23 # R, 128 bytes function(x,m,n)t(round(which((1:m*det(x))%%m<1.5)*outer(1:n,1:n,Vectorize(function(a,b)det(x[-a,-b,drop=F])*(-1)^(a+b))))%%m)  Try it online! A function taking three arguments, x = the matrix, m = the modulus and n the number of rows of x. Returns a matrix. Uses the same method as my Jelly answer. # Jelly, (21?) 22 bytes The trailing Ṫ may be unnecessary - it is if there is only ever a single inverse of M with elements modulo m. Ḷṗ⁹L²¤ṁ€⁹æ×%³L⁼þ$ƑɗƇṪ


A full program printing the result.

Try it online! - Too slow for any of the given test cases (the 35 case took ~20 minutes locally).

11 bytes (but floating point output):

Using Bubler's observation (go upvote!) that raising the determinant to Euler's totient is enough to remove the determinant's denominators:

æ*-×ÆḊ*ÆṪ}ɗ


However, unlike in J, the inversion of $$\M\$$ in Jelly gives floats so we no longer get an integer matrix as output.

Try it online!

• @RGS - I guess the 11 would only be acceptable if an integer output was given, even though this is not explicitly stated in the question? Mar 10 '20 at 20:52

# WIN+APL, 114 bytes

Prompts for matrix followed by modulus.

m←r←⎕⋄z←r[1;1]⋄⍎∊(¯1+1↑⍴r)⍴⊂'z←z×1 1↑r←(1 1↓r)-((1↓r[;1])∘.×1↓r[1;])÷r[1;1]⋄'⋄⌊.5+n|((1=n|z×⍳n)/⍳n←⎕)×(z←⌊.5+z)×⌹m


Try it online! Courtesy of Dyalog Classic

# Magma, 34 bytes

func<m,M|Matrix(Integers(m),M)^-1>


No TIO for magma, though you can try it on http://magma.maths.usyd.edu.au/calc/

# Java 8, 270 261 bytes

M->m->{int l=M.length,R[][]=new int[l][l],T[][]=new int[l][l],d=0,s=l,r,c,k;for(;d!=1|s!=0;){for(r=l*l;r-->0;R[r/l][r%l]=d*=Math.random())d=m;for(d=1,s=r=l;r-->0;d*=T[r][r]%m)for(c=l;c-->0;s-=T[r][c]%m)for(T[r][c]=k=0;k<l;)T[r][c]+=M[r][k]*R[k++][c];}return R;}


-9 bytes thanks to @ceilingcat.

Keeps trying random matrices (including duplicates) until it find the correct one, so times out for most test cases. I tried adding a cache so it tries random matrices without duplicates, but then it still times out for the same test cases.

Try it online (only contains the test cases m=35; M=[[24,14],[48,45]] and m=5; M=[[15,13],[21,13]]).

Explanation:

M->m->{                    // Method with int-matrix & int parameters and int-matrix return
int l=M.length,          //  Dimension of the input-matrix
R[][]=new int[l][l], //  Result-matrix of that same size
T[][]=new int[l][l], //  Temp-matrix of that same size
d=0,                 //  Flag for the diagonal
s=l,                 //  Flag for the decreasing sum
r,c,k;               //  Index integers
for(;d!=1                //  Continue looping as long as the diagonal flag isn't 1 yet
|s!=0;){            //  nor the decreasing sum flag isn't 0 yet:
for(r=l*l;r-->0;       //   Loop over all cells:
R[r/l][r%l]=         //     Set the current cell in matrix R:
d*=Math.random())d=m;
//      To a random value in the range [0,m)
for(d=1,               //   Reset the diagonal flag to 1
s=r=l;             //   Reset the decreasing sum flag to l
r-->0              //   Loop over the rows:
;                  //     After every iteration:
d*=               //      Multiply the diagonal flag by:
T[r][r]        //       The value in the r,r'th cell of matrix T
%m)     //       Modulo the input m
for(c=l;c-->0        //    Inner loop over the columns:
;                //      After every iteration:
s-=             //       Decrease the decreasing sum flag by:
T[r][c]      //        The value in the r,c'th cell of matrix T
%m)   //        Modulo the input m
for(T[r][c]=k=0;   //     Reset the r,c'th cell of matrix T to 0
k<l;)          //     Inner loop k in the range [0, length):
T[r][c]+=        //      Increase the r,c'th cell of matrix T by:
M[r][k]        //       The r,k'th cell of matrix M
*R[k++][c];}   //       Multiplied by the k,c'th cell of matrix R
return R;}               //  And if the loops are done: return matrix R as result


# R, 97 83 bytes

function(M,m,d){while(any(M%*%(x=matrix(T%/%m^(1:d^2-1),d))%%m-diag(d)))T=T+1;x%%m}


Try it online!

Pretty slow. Takes the dimension of the matrix as input. The previous version using a for loop is a bit faster.

Thanks to Robin Ryder for -14 bytes.

### Explanation:

We iterate over every number between $$\1\$$ and $$\m^{d^2}\$$, converting each to its base-$$\m\$$ digits (with leading zeros), reshaping those digits into a matrix of the appropriate size, and testing to see if it's the inverse of $$\M\$$ modulo $$\m\$$.

I wanted to attempt the whole series in SNOBOL but I'm not sure I will be able to implement matrix multiplication in SNOBOL in time for it to be a valid submission...

• 83 bytes with a while loop. Mar 10 '20 at 16:33

# Python 3 + SymPy, 33 bytes

from sympy import*
Matrix.inv_mod
`

Try it online!

SymPy's Matrix class has a method for modular inverse.