Ask Your Question
2

Solving systems of congruences with Gaussian integers

asked 2018-11-15 03:57:43 -0600

Jsevillamol gravatar image

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?

edit retag flag offensive close merge delete

Comments

1

Note that CRT is a theorem for finitely generated abelian groups. That the additive group Z[i] also has a (euclidean) ring structure is not required for CRT. All you need is a list of finite index subgroups that are "pairwise coprime" (that any two of then generate the entire additive group). In particular, it works with any ring of integer and a list of coprime ideals. It doesn't need the ideals to be principal. See idealchinesein pari/gp: https://pari.math.u-bordeaux.fr/docht...

nbruin gravatar imagenbruin ( 2018-11-15 18:03:41 -0600 )edit

1 answer

Sort by ยป oldest newest most voted
2

answered 2018-11-15 13:41:44 -0600

rburing gravatar image

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
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-11-15 03:57:43 -0600

Seen: 213 times

Last updated: Nov 15 '18