(RGS 4/5) Inverting matrices modulo m

Question

Task

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 code-golf 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 code-golf challenge, so enjoy!

Answer 1

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

Answer 2

J, ¹⁸ 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)

Answer 3

Wolfram Language (Mathematica), 23 bytes

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

Inverse[#2,Modulus->#]&

Try it online!

Answer 4

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
    0                        //   start with s = 0
  )                          // 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()

Answer 5

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

Answer 6

SageMath, ⁴⁸ 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.

Answer 7

Sledgehammer, 6 bytes

⠑⡿⡆⠱⣁⣭

Decompresses into this Wolfram Language function:

Inverse[#2, Modulus -> #1]

Try it online!

Answer 8

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.

Answer 9

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)

Answer 10

R, 128 bytes

function(x,m,n)t(round(which((1:m*det(x))%%m<1.5)[1]*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.

Answer 11

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!

Answer 12

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

Answer 13

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/

Answer 14

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

Answer 15

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

Answer 16

Python 3 + SymPy, 33 bytes

from sympy import*
Matrix.inv_mod

Try it online!

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