I am writing my own function to solve systems of congruences in arbitrary euclidean domains, but I can't quite make it work for gaussians.

```
# Extended Euclides Algorithm
def extended_euclides(a,b,quo=lambda a,b:a//b):
r0 = a; r1 = b
s0 = 1; s1 = 0
t0 = 0; t1 = 1
while r1 != 0:
q = quo(r0, r1)
r0, r1 = r1, r0 - q * r1
s0, s1 = s1, s0 - q * s1
t0, t1 = t1, t0 - q * t1
return r0, s0, t0
# Chinese Remainder Theorem
def chinese_remainder(remainders, modules, quo=lambda a,b:a//b):
"""
INPUT:
A list of remainders v_i and pairwise coprime modules m_i
representing a system of congruences {x == v_i mod m_i},
where v_i, m_i belong to an euclidean domain ED
OUTPUT:
A solution of the system of congruences in the domain ED
"""
m = reduce(lambda x,y:x*y, modules)
c = 0
for v_i, m_i in zip(remainders, modules):
m_div_m_i = quo(m, m_i)
_, s_i, _ = extended_euclides(m_div_m_i, m_i, quo)
# c_i = (v_i * s_i) % m_i
# we need to do it convolutedly bcs gaussian integer quotient is not implemented in SageMath
a = v_i * s_i; b = m_i
c_i = a - quo(a,b)*b
c += c_i * m_div_m_i
return c
# EXAMPLE ON Z
remainders = [ZZ(2),ZZ(7)]
modules = [ZZ(11), ZZ(13)]
print("System to be solved:")
for r,m in zip(remainders, modules):
print("x == {} (mod {})".format(r,m))
solution = chinese_remainder(remainders, modules)
print("A solution of the system is {}".format(solution))
#check that the solution is correct
for r, m in zip(remainders, modules):
assert(solution % m == r)
# EXAMPLE ON Z[i]
ZI = QuadraticField(-1, 'I').ring_of_integers()
gaussian_quo = lambda a,b : ZI(real(a/b).round() + I*imag(a/b).round())
gaussian_rem=lambda a,b: a - gaussian_quo(a,b)*b
remainders = [ZI(2+3*I),ZI(7+5*I)]
modules = [ZI(11), ZI(13-I)]
print("System to be solved:")
for r,m in zip(remainders, modules):
print("x == {} (mod {})".format(r,m))
solution = chinese_remainder(remainders, modules, gaussian_quo)
print("A solution of the system is {}".format(solution))
#check that the solution is correct
for r, m in zip(remainders, modules):
print(gaussian_rem(solution, m) / r)
assert(gaussian_rem(solution, m) == r)
```

This produces the following ouptut:

```
System to be solved:
x == 2 (mod 11)
x == 7 (mod 13)
A solution of the system is 46
System to be solved:
x == 3*I + 2 (mod 11)
x == 5*I + 7 (mod -I + 13)
A solution of the system is -41*I - 75
1
36/37*I - 6/37
Error in lines 31-33
Traceback (most recent call last):
File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute
flags=compile_flags) in namespace, locals
File "", line 3, in <module>
AssertionError
```

So in the integer example we find a correct solution, but in the Gaussian example we find a solution that satisfies the first incongruence relation, but not the second.

I have checked with wolfram alpha, and if I did not mess up it seems like the problem is in the Chinese Remainder Theorem algorithm rather than on the assertion tests.

What is going on?