Solving matrices in a finite field

[Andrew Mizener]

For my project, I need to solve for a matrix X given matrices Y and K. (X*Y=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?