Easiest way to perform modular matrix inversion with Python?

Question

I'd like to take the modular inverse of a matrix like [[1,2],[3,4]] mod 7 in Python. I've looked at numpy (which does matrix inversion but not modular matrix inversion) and I saw a few number theory packages online, but nothing that seems to do this relatively common procedure (at least, it seems relatively common to me).

By the way, the inverse of the above matrix is [[5,1],[5,3]] (mod 7). I'd like Python to do it for me though.

Answer 1

A hackish trick which works when rounding errors aren't an issue:

• find the regular inverse (may have non-integer entries), and the determinant (an integer), both implemented in numpy
• multiply the inverse by the determinant, and round to integers (hacky)
• now multiply everything by the determinant's multiplicative inverse (modulo your modulus, code below)
• do entrywise mod by your modulus

A less hackish way is to actually implement gaussian elimination. Here's my code using Gaussian elimination, which I wrote for my own purposes (rounding errors were an issue for me). q is the modulus, which is not necessarily prime.

def generalizedEuclidianAlgorithm(a, b):
    if b > a:
        return generalizedEuclidianAlgorithm(b,a);
    elif b == 0:
        return (1, 0);
    else:
        (x, y) = generalizedEuclidianAlgorithm(b, a % b);
        return (y, x - (a / b) * y)

def inversemodp(a, p):
    a = a % p
    if (a == 0):
        print "a is 0 mod p"
        return None
    if a > 1 and p % a == 0:
        return None
    (x,y) = generalizedEuclidianAlgorithm(p, a % p);
    inv = y % p
    assert (inv * a) % p == 1
    return inv

def identitymatrix(n):
    return [[long(x == y) for x in range(0, n)] for y in range(0, n)]

def inversematrix(matrix, q):
    n = len(matrix)
    A = np.matrix([[ matrix[j, i] for i in range(0,n)] for j in range(0, n)], dtype = long)
    Ainv = np.matrix(identitymatrix(n), dtype = long)
    for i in range(0, n):
        factor = inversemodp(A[i,i], q)
        if factor is None:
             raise ValueError("TODO: deal with this case")
        A[i] = A[i] * factor % q
        Ainv[i] = Ainv[i] * factor % q
        for j in range(0, n):
            if (i != j):
                factor = A[j, i]
                A[j] = (A[j] - factor * A[i]) % q
                Ainv[j] = (Ainv[j] - factor * Ainv[i]) % q
    return Ainv

EDIT: as commenters point out, there are some cases this algorithm fails. It's slightly nontrivial to fix, and I don't have time nowadays. Back then it worked for random matrices in my case (the moduli were products of large primes). Basically, the first non-zero entry might not be relatively prime to the modulus. The prime case is easy since you can search for a different row and swap. In the non-prime case, I think it could be that all leading entries aren't relatively prime so you have to combine them