Revision history [back]

The code and output is correct. The remainder after dividing -41*I - 75 by -I + 13 is 6*I - 6, which is congruent to 5*I + 7 modulo -I + 13. Your final comparison

gaussian_rem(solution, m) == r

is not the right thing to test, because it implicitly assumes that the right-hand sides of your system of equations were already in normal form. Instead, you want to test

gaussian_rem(solution, m) == gaussian_rem(r, m)

or more efficiently

gaussian_rem(solution - r, m) == 0