[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](/questions/tagged/code-golf "show questions tagged '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](https://codegolf.meta.stackexchange.com/q/18545/75323). 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](https://codegolf.meta.stackexchange.com/q/18545/75323))
Also, as per section 4 of the rules in the [linked meta post](https://codegolf.meta.stackexchange.com/q/18545/75323), the "restricted languages" for this third challenge are *only* [Jelly](https://codegolf.stackexchange.com/a/200519/75323), [V (vim)](https://codegolf.stackexchange.com/a/200525/75323) and [05AB1E](https://codegolf.stackexchange.com/a/200491/75323) 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](/questions/tagged/code-golf "show questions tagged 'code-golf'") challenge, so enjoy!
  
[Answer]
# [R](https://www.r-project.org/), 68 bytes
```
function(M,m,n,A=M){while(any(A%*%M%%m!=diag(n)))A[]=rpois(n^2,9)
A}
```
[Try it online!](https://tio.run/##K/qfZvs/rTQvuSQzP0/DVydXJ0/H0dZXs7o8IzMnVSMxr1LDUVVL1VdVNVfRNiUzMV0jT1NT0zE61raoID@zWCMvzkjHUpPLsfZ/mkZuYklRZoVGsoaRmY6CIRAbW@gomBhq6hhpAmlTHQUjzf8A "R – 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]
# [J](http://jsoftware.com/), ~~18~~ 16 bytes
```
(]%1+.]^5 p:[)%.
```
[Try it online!](https://tio.run/##dVBNT8JAED3TXzEx9GOLLu62Fd0Es9HEk/HgtWmNgRLrRyCAZo2Vv44zS4GC6WGy2Tdv3nszr@uH6bIAfzYvvsrp5wKW87JY@E4MyjdVGY5DASVfmSoouWHheMhD3j8NxI/0jGLedTnkLv/2nbRyuc629LQquU6JnSrk68aAdnkrWViu/uWyUtooyzzrAw@RbCd23/xJ9BKYqRQZ400vQIG6y/btIHN7VjTLLcBQFLFUWbQJBmnNVDlLvJliK0K3YJBbzIICb@NyN/cCmv5gOrAchqv5Vkf0eFMjI6AWWNUBhIecZiICjkLSgE7zFU5nBG3@lH6H6cxhzgn3JzBU4MMpnIPCOuNw@3h/t95Y71XXzPl4XhJXguzKCwMCK7o0EAvjFKOXKcQJTABJnU7n4YbDvFh8vi93rYpa0KMz75nEG70Uo7cD8SsDybkBerfa0VWrNraa2jvmP@0Ioq5IUFNgcNLHEhIXwDeK8I89ud2l3S8@9Itb/WKIu3gf3EKiixjgG2@cqQTiERZdkIwxQFLnkjZPvd2gffHB4eKD4yDrPw "J – 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](http://jsoftware.com/), 18 bytes
```
%.@]*-/ .*@]^5 p:[
```
[Try it online!](https://tio.run/##dZBbT8JAEIWf6a@YGKC0wMBeamUTTKOJT4YHXxs1BkpYL4FAMWus/PU6W2gBTR8mk93z7Tkz@5pPlmkC7mqdfOrldgPpWicb15GgXJNpf@Yz0LgzWUej8fzZGH0c9Drsm7eN8trXeowt/HKdOGth9FjicaYxii0dK@KjkwdRC2thVrDRD/JMRUYVZH8A6BNcvKiOT8@sG8BKxUTM9lqHDA6qd5Q95wLdOYwVuNCDISiqPsLtw/1dXjwrDQs89xzn4yW1OAfe5JcGGJW4MiCZcZLpYgkygDkQ1Gg0JjcI62SzfU8rKbMSdK3pkbTcdJFM387MRwaCoQHbS28xqvUm6dS7Iv95CxBNFpAno8GtPxXjtAB1IehMGi93qc@T53myNk@CbNL/0BacUlhIXe6TbTG6F1T2B20wDRAc5uLFPIftwvrFw/PFw7@D5L8 "J – 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]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 23 bytes
¯\\_(ツ)\_/¯ the answer was in the documentation of [Modulus](http://reference.wolfram.com/language/ref/Modulus.html)
```
Inverse[#2,Modulus->#]&
```
[Try it online!](https://tio.run/##y00syUjNTSzJTE78n26r8N8zryy1qDg1WtlIxzc/pTSntFjXTjlW7X9AUWZeiUN6tLG5TnW1hY6ljpGRjqF5rU61kYmOsQEIGVoCecaWOhY6JqY6RsZAjilIGKjO2Li2NvY/AA "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# 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!](https://tio.run/##pZNNc5swEIbv/hW6tKyKIAhBHLsDuaRHTjlSTcLY2CYY8MiUSWby391dyWSmp04bDxZarfQ@@yFeqqk6b0xzGoN@2NaXXXbpsrzAJ@yqE4ARbzzLjTOexCsaAPusxZfvt98eoOBfOr/DIZD3e2j5@tV/@yrvg3bd8udn3HCAAo9ZAZPlB6vI50NcHDKoRI@yVbhrjmNt4CnL@yDg4iHDMPzi/b0ITb39takBzqKctGhw99mHhjDTeuJXSPMHJEIIDZfv5YKVLEkFK8v4VjB5q3Gq7gRLpNYMjX/93dyQVIwCiqTilWBqqTVh1MpicEyj2fefFIeRFPGKpBKkpFeKSwZXZGJdlEz6iWRQbzkzlGOkihgJjrHNROKeu89kks7RqiXRLCSx5SJpKRHtShahRcW1lqIA0B//FW0hdjtWI3adIWG0otmgQslrghRGWeICxkDNlLYCVFQVub@rO7XU1ndWTZ2XzlCxmJ5bReIoSpcs/ui@IkvRLaAi4kJkiVbJypAPoTFek4VehLvB/Kg2B4CyE6zQnGU52wz9eTjW4XHYww46jh@du@rkdN/mRNMpHIfH0TT9Hnh4qraPY2VGiDnzmcc8Hr4MTQ@ex@fZz96zPnrzy28 "JavaScript (Node.js) – 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]
# [Jelly](https://github.com/DennisMitchell/jelly), 25 bytes
```
ÆḊ×Ɱ⁹%ỊTḢ×ZÆḊ-Ƥ$-ƤNÐe⁺€Zʋ
```
[Try it online!](https://tio.run/##y0rNyan8//9w28MdXYenP9q47lHjTtWHu7tCHu5YdHh6FFhc99gSFSD2Ozwh9VHjrkdNa6JOdf8/vByk8v//6GhDUx0FQ0MdBWODWB2FaCMDIM9IR8EEzDM21lEwAsobxcb@N7EEAA "Jelly – 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]
# [SageMath](https://doc.sagemath.org/html/en/developer/index.html), ~~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`.
[Answer]
# [Sledgehammer](https://github.com/tkwa/Sledgehammer), 6 bytes
```
⠑⡿⡆⠱⣁⣭
```
Decompresses into this Wolfram Language function: 
```
Inverse[#2, Modulus -> #1]
```
[Try it online!](https://tio.run/##y00syUjNTSzJTE78n2ar8N8zryy1qDg1WtlIR8E3P6U0p7RYQddOQdkwVu1/QFFmXolDWrSxuU51tYWOpY6RkY6hea1OtZGJjrEBCBlaAnnGljoWOiamOkbGQI4pSBiozti4tjb2PwA "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# [Charcoal](https://github.com/somebody1234/Charcoal), 41 bytes
```
ＦＥＸθ×ηη⪪Ｅ×ηη÷ιＸθλη¿⬤ι⬤ζ⁼⁼λν﹪ΣＥμ×ξ§§ιπλθＩι
```
[Try it online!](https://tio.run/##TY09D4IwEIb/yo3X5BwwOjkZdXAgIcGNMDRQQpOD8lGQ@Odri2Kcnsvd@9xb1HIojGTnKjMAxrLDxDzVgD3BQzdqxJqgFoIg7VjbNfC3J7i39qpnXSrUBD@VRVC8J0BXgGfmcA54Edz6SfKIXzBB67OxKSc2mE7N2tFs9YvX7L0t1YIb/adOfDo8@gBIBt1avMjRohbi5Fx2JNgTZFnkmfvhQBDlee52M78B "Charcoal – 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]
# [MATL](https://github.com/lmendo/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!](https://tio.run/##y00syfn/3yozLypOUUnN3SE6NjVSqzymJKoyotLW3kHv/39TrmhDUwVDY2sFI0MgFQsA "MATL – 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]
# [R](https://www.r-project.org/), 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!](https://tio.run/##VY7NasMwDMfvewpfClKmjDq222Ys1z3CLqWD5qvNIUkxCS17@Uyyx2gOAlv6/T/8cimWdh6qqRsHeFBPA07gx3mo4X7tqiuAfu@TupnggbjZ9B/6zeFRn5JxnhrPx4FkvppqGn3308C/2ZlKDLpjeqa0pNqPt@LzhAmkGr/h/FoiBktc2uK5Al5iET/eJfSlBV923KeCbEdK75BUBeZAympEsm5F5KTcNhDyFMDkK8Dy1gbAioUTYmVhDaksWmhH6sCAM8@AbLUmZf5ytvzN2Cp@jciZyKTbKprjuFLGqN5HpRWTODqPaiZCLQZN2Lh4F5kxUnaPyy8 "R – 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](https://codegolf.stackexchange.com/a/200691/42248). 
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), (21?) 22 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
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!](https://tio.run/##y0rNyan8///hjm0Pd05/1LjT59CmQ0se7mx81LQGyDu87PB01UObfR417jm8L0Hl2MST04@1P9y56v///6b/o6MNTXUUDI1jdRSijQzBrFgA "Jelly – 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](https://codegolf.stackexchange.com/a/200834/53748) (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!](https://tio.run/##y0rNyan8///wMi3dw9MPtz3c0aUFJHeuqj05/f///9HRhqY6CoaGOgrGBrE6CtFGBkCekY6CCZhnbKyjYASUN4qN/W9iCQA "Jelly – Try It Online")
[Answer]
# 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](https://tio.run/##PU7BSsNAEL3vV8ytu5RUZ3fTNhT/oP2C0kOoRIQ0Snoy9CTShrQpikj9Cw8K4tX8yfxInElB2Bnee/vmzcT3aXD9EKd3N8Eyjdfr22VLx7fZlLbPTlG5S9oVw5yLZdo/FcLmOMEFE6qPVFb69wP7SNsXqr9yw42qx574iuaEIB8yrjv4mptA6w7MOcNQ@T5gV8dxsjDN9394TxYcqkHYzzY8cpVtOI/qz8xcSD9fZJqTllVnY8GUDj@rlg9vE2XB8jV2CDgENwaPyocqUR48y2OIwFrAEVgP7lIeRuAiYGMI1kEoEjucU270Bw "APL (Dyalog Classic) – Try It Online")
[Answer]
# [Magma](http://magma.maths.usyd.edu.au/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]
# 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](https://tio.run/##hZJRa9swEMff8ynUQMGKZS1uExhVZdjLoA9mkObN@EGzndSxJBtJ7gieP3t2dtWGQdnA4PPd7//X6c4n8SqiU9lcCimsRamo9bBAyDrh6gKdoEp7V0t66HXh6lbT7z54rLXL8iwn/4SetKuOlSHI00mCDohf0ihRUTJAEkmeUlnpo3shu4nguvo10zKHh@w/yZV8TSyXxJCCNOzQmoCVNzz@bW/4muFhShguV5KZKErWbJeZLyA0tzLn5Yqnwr1QI3TZqgDjkqvZoeQxeILMiwDcZwZU@a3CE1BAqZhLNnorFe8l/8UbOL55lAz7RMjTKWjy1S5rwnBKsdFUrjca7dh4YQuYdNf/lDBpP/DXti6RgiUEz87U@pjlSOBpIQi5yrpgSxB6HwYMZhjiLYnvRzLcxdN7xOzK3gP8F3u3IfEG2M1Xstl6dlxctz0fPkunxaiPpaHUt/B8tq5StO0d7aA7J3WwfNJd7@wDUnwZqnDJUArB9Zf4Zow4W1pWVbdv3@4UpNi3@Znfj96B4QP6n8mBiq6TZzDzgcLY32i8/AE) (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]
# [R](https://www.r-project.org/), ~~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!](https://tio.run/##TczRCsIwDAXQd79iIIVUb5Fu3RB1n7C3Pg/G5rRgJ4yJFfHbawYKgxDCzUnG2CcnlcT@MbSTuw9UwaOT7@fV3c7UDC@qxEZQKH0zjS6QFTvha9KHrk6VlkylEF51rrkQz9KWdquPgbNPXPf0u2pJI0MOI5FyFdxXyy0nukC2h@GXKWAlTD6jpTHQBoZN/jfZbOIX "R – Try It Online")
Pretty slow. Takes the `d`imension 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]
# [Python 3](https://docs.python.org/3/) + [SymPy](https://docs.sympy.org/latest/index.html), 33 bytes
```
from sympy import*
Matrix.inv_mod
```
[Try it online!](https://tio.run/##XVLLaoRAELz7FQ256DIJzrSzroHkD/IFMoRAkHjwgWtC/HrTNY2b1YNKd3VVT9U4LvPX0PO6NtPQ0XXpxoXabhym@ZS8fcxT@/vU9j/v3fC5Ni/7RpI0w0RtP37PV/lQnRClOpLWtTsbsudgqOaLocKGkMnHy@uBHl9JBpzUjAFXGeIyhD1fmj7fYKVz9U@30K@AFwL4I1uatogotnul3233hsqNzAdywYZcXG1l7BK5nndc1eYSKnsyKNaKqp49lwpGY8UQFtxpGHduIiIndRoINKTKtwIm7PGY0pNkEKONVuCZc300GRY82t@EvaLgwDUyKbNbolgijnBz7pY8o2LcAPKQRh7XRrmoBUw2O7miJDzTOLX9nDbpSX@MLFv/AA "Python 3 – Try It Online")
SymPy's Matrix class has a [method for modular inverse](https://docs.sympy.org/latest/modules/matrices/matrices.html?highlight=inv_mod#sympy.matrices.matrices.MatrixBase.inv_mod).
]