How to solve a congruence system in python?

Question

for the problem Ax ≡ B (MOD C) I did this, and it was okay:

 def congru(a,b,c):
    for i in range(0,c):
       if ((a*i - b)%c)== 0 :
          print(i)

Now I have to solve a system of equations, where A = ( 5x + 7y) and A= (6x + 2y), and B= 4 and B = 12 , respectively, and C is 26.

In other words: ( 5x + 7y)≡ 4 (mod 26) (6x + 2y)≡ 12 (mod 26)

How do I do that?

Thanks.

Answer 1

For the aX ≡ b (mod m) linear congruence, here is a more powerful Python solution, based on Euler's theorem, which will work well even for very large numbers:

def linear_congruence(a, b, m):
    if b == 0:
        return 0

    if a < 0:
        a = -a
        b = -b

    b %= m
    while a > m:
        a -= m

    return (m * linear_congruence(m, -b, a) + b) // a

>>> linear_congruence(80484954784936, 69992716484293, 119315717514047)
>>> 45347150615590

For the X, Y system: multiply 1st equation by 2 and the 2nd one by -7, add them together and solve for X. Then substitute X and solve for Y.

Answer 2

This is a very fast linear congruence solver that can solve a 4096 byte number in a second. Recursion versions meet their max depth at this length:

#included if you use withstats=True, not needed otherwise
def ltrailing(N):
    return len(str(bin(N))) - len(str(bin(N)).rstrip('0'))

#included if you use withstats=True, not needed otherwise
def _testx(n, b, withstats=False):
    s = ltrailing(n - 1)
    t = n >> s
    if b == 1 or b == n - 1:
        return True
    else:
        for j in range(0, s):
            b = pow_mod(b, 2, n, withstats)
            if b == n - 1:
                return True
            if b == 1:
                return False
            if withstats == True:
               print(f"{n}, {b}, {s}, {t}, {j}")
    if withstats == True:
        print(f"{n}, {b}, {s}, {t}")
    return False

def fastlinearcongruence(powx, divmodx, N, withstats=False):  
 x, y, z = egcditer(powx, N, withstats)  
 answer = (y*divmodx)%N  
 if withstats == True:  
    print(f"answer = {answer}, mrbt = {a}, mrst = {b},  mranswer = {_testx(N, answer)}")  
 if x > 1:  
    powx//=x  
    divmodx//=x  
    N//=x  
    if withstats == True:  
      print(f"powx = {powx}, divmodx = {divmodx}, N = {N}")  
    x, y, z = egcditer(powx, N, withstats)  
    if withstats == True:  
      print(f"x = {x}, y = {y}, z = {z}")  
    answer = (y*divmodx)%N  
    if withstats == True:  
       print(f"answer = {answer}, mrbt = {a}, mrst = {b},  mranswer = {_testx(N, answer)}")  
 answer = (y*divmodx)%N  
 if withstats == True:  
    print(f"answer = {answer}, mrbt = {a}, mrst = {b},  mranswer = {_testx(N, answer)}")  
 return answer 
 
def egcditer(a, b, withstats=False):  
  s = 0  
  r = b  
  old_s = 1  
  old_r = a  
  quotient = 0 
  if withstats == True:  
      print(f"quotient = {quotient}, old_r = {old_r}, r = {r}, old_s = {old_s}, s = {s}")  
  while r!= 0:  
    quotient = old_r // r  
    old_r, r = r, old_r - quotient * r  
    old_s, s = s, old_s - quotient * s  
    if withstats == True:  
      print(f"quotient = {quotient}, old_r = {old_r}, r = {r}, old_s = {old_s}, s = {s}")  
  if b != 0:  
    bezout_t = quotient = (old_r - old_s * a) // b  
    if withstats == True:  
      print(f"bezout_t = {bezout_t}")  
  else:  
    bezout_t = 0  
  if withstats == True:  
    print("Bézout coefficients:", (old_s, bezout_t))  
    print("greatest common divisor:", old_r)  
  return old_r, old_s, bezout_t

To use:

In [2703]: fastlinearcongruence(63,1,1009)                                                                                                                                                   
Out[2703]: 993


In [2705]: fastlinearcongruence(993,1,1009)                                                                                                                                                  
Out[2705]: 63


Also this is 100x faster than pow when performing this pow: 

pow(1009, 2**bit_length()-1, 2**bit_length())

where the answer is 273.

This is equivalent to the much faster:


In [2713]: fastlinearcongruence(1009,1,1024)                                                                                                                                              
Out[2713]: 273