Sympy: Solving Matrices in a finite field

Question

For my project, I need to solve for a matrix X given matrices Y and K. (XY=K) The elements of each matrix must be integers modulo a random 256-bit prime. My first attempt at solving this problem used SymPy's mod_inv(n) function. The problem with this is that I'm running out of memory with matrices of around size 30. My next thought was to perform matrix factorization, as that might be less heavy on memory. However, SymPy seems to contain no solver that can find matrices modulo a number. Any workarounds or self-made code I could use?

Raspings

Answer 1

sympy's Matrix class supports modular inverses. Here's an example modulo 5:

from sympy import Matrix, pprint

A = Matrix([
    [5,6],
    [7,9]
])

#Find inverse of A modulo 26
A_inv = A.inv_mod(5)
pprint(A_inv)

#Prints the inverse of A modulo 5:
#[3  3]
#[    ]
#[1  0]

The rref method for finding row-reduced echelon form supports a keyword iszerofunction that indicates what entries within a matrix should be treated as zero. I believe the intended use is for numerical stability (treat small numbers as zero) although I am unsure. I have used it for modular reduction.

Here's an example modulo 5:

from sympy import Matrix, Rational, mod_inverse, pprint

B = Matrix([
        [2,2,3,2,2],
        [2,3,1,1,4],
        [0,0,0,1,0],
        [4,1,2,2,3]
])

#Find row-reduced echolon form of B modulo 5:
B_rref = B.rref(iszerofunc=lambda x: x % 5==0)

pprint(B_rref)

# Returns row-reduced echelon form of B modulo 5, along with pivot columns:
# ([1  0  7/2  0  -1], [0, 1, 3])
#  [                ]
#  [0  1  -2   0  2 ]
#  [                ]
#  [0  0   0   1  0 ]
#  [                ]
#  [0  0  -10  0  5 ]  

That's sort of correct, except that the matrix returned by rref[0] still has 5's in it and fractions. Deal with this by taking the mod and interpreting fractions as modular inverses:

def mod(x,modulus):
    numer, denom = x.as_numer_denom()
    return numer*mod_inverse(denom,modulus) % modulus

pprint(B_rref[0].applyfunc(lambda x: mod(x,5)))

#returns
#[1  0  1  0  4]
#[             ]
#[0  1  3  0  2]
#[             ]
#[0  0  0  1  0]
#[             ]
#[0  0  0  0  0]

Culosio

NB: This function does not always work. An example is Matrix([[4,3,1,3],[2,4,1,3]]) in Z_5. In this case, the regular iszerofunc call using lambda x: x % 5 == 0 gives a matrix with denominators including a 5. Since in Z_5 there is no inverse for 5, the program will exit. – Trumpery

@Trumpery In Sympy 1.6, this example now works! (It fails as you describe in Sympy 1.0). I never looked into why this example failed in 1.0, so I cannot say for certain that all examples will now work, but happily the method seems more reliable now. (I'm also not sure what version this started working.) – Culosio