Solving systems of congruences with Gaussian integers

savecancel

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?