I'm looking for a better algorithm than one I found on stackoverflow to handle 4096 byte numbers, i'm hitting a maximum recursion depth.
Code from stackoverlow post, taken from this Stack Overflow post:
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
This works fine for smaller numbers, for example:
In [167]: pow_mod(8261, 63, 4033)
63 1 8261 4033
31 195 1728 4033
15 2221 1564 4033
7 1231 2098 4033
3 1518 1601 4033
1 2452 2246 4033
0 2147 3266 4033
Out[167]: 2147
And the linear congruence works:
linear_congruence(8261, 3266, 4033):
2147
But i hit maximum recursion depth with larger numbers. Is there a better algorithm or non recursive algorithm of the linear_congruence algorithm i provided?
Based on Eric Postpischil's remark, i wrote the pseudocode from the wikipedia entry and created a very fast linear congruence algorithm utilizing the method from here: http://gauss.math.luc.edu/greicius/Math201/Fall2012/Lectures/linear-congruences.article.pdf .
This works nicely on pows with a powers of 2-1, to get the answer. I'm looking into how offsetting from this changes the answer and hope to incorporate it to work for those answers as well, but for now, i have what i need since i'm working with powers of 2 -1 for y in pow(x, y, z):
def fastlinearcongruencex(powx, divmodx, N, withstats=False):
x, y, z = egcditerx(powx, N, withstats)
if x > 1:
powx//=x
divmodx//=x
N//=x
if withstats == True:
print(f"powx = {powx}, divmodx = {divmodx}, N = {N}")
x, y, z = egcditerx(powx, N)
if withstats == True:
print(f"x = {x}, y = {y}, z = {z}")
answer = (y*divmodx)%N
if withstats == True:
print(f"answer = {answer}")
return answer
def egcditerx(a, b, withstats=False):
s = 0
r = b
old_s = 1
old_r = a
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
It even works instantaneously on 4096 byte numbers which is great:
In [19036]: rpowxxxwithbitlength(1009,offset=0, withstats=True, withx=True, withbl=True)
63 1 272 1009
31 272 327 1009
15 152 984 1009
7 236 625 1009
3 186 142 1009
1 178 993 1009
0 179 256 1009
Out[19036]: (179, 256, True, 272)
In [19037]: fastlinearcongruencex(272,256,1009)
Out[19037]: 179
Thank you Eric for pointing out what this was, i wrote an extremely fast linear congruence algorithm utilizing egcd and the procedure from the pdf above. If any stackoverflowers need a fast algorithm, please point them to this one. I also learned that congruence is always maintained when the pow(x,y,z) has a y off of the powers of 2-1. I will look into this further to see if an offset change exists to keep the answers intact and will follow up in the future if found.
If you have Python 3.8 or later, you can do everything you need to with a very small number of lines of code.
First some mathematics: I'm assuming that you want to solve ax = b (mod m) for an integer x, given integers a, b and m. I'm also assuming that m is positive.
The first thing you need to compute is the greatest common divisor g of a and m. There are two cases:
if b is not a multiple of g, then the congruence has no solutions (if ax + my = b for some integers x and y, then any common divisor of a and m must also be a divisor of b)
if b is a multiple of g, then the congruence is exactly equivalent to (a/g)x = (b/g) (mod (m/g)). Now a/g and m/g are relatively prime, so we can compute an inverse to a/g modulo m/g. Multiplying that inverse by b/g gives a solution, and the general solution can be obtained by adding an arbitrary multiple of m/g to that solution.
Python's math module has had a gcd function since Python 3.5, and the built-in pow function can be used to compute modular inverses since Python 3.8.
Putting it all together, here's some code. First a function that finds the general solution, or raises an exception if no solution exists. If it succeeds, it returns two integers. The first gives a particular solution; the second gives the modulus that provides the general solution.
def solve_linear_congruence(a, b, m):
""" Describe all solutions to ax = b (mod m), or raise ValueError. """
g = math.gcd(a, m)
if b % g:
raise ValueError("No solutions")
a, b, m = a//g, b//g, m//g
return pow(a, -1, m) * b % m, m
And then some driver code, to demonstrate how to use the above.
def print_solutions(a, b, m):
print(f"Solving the congruence: {a}x = {b} (mod {m})")
try:
x, mx = solve_linear_congruence(a, b, m)
except ValueError:
print("No solutions")
else:
print(f"Particular solution: x = {x}")
print(f"General solution: x = {x} (mod {mx})")
Example use:
>>> print_solutions(272, 256, 1009)
Solving the congruence: 272x = 256 (mod 1009)
Particular solution: x = 179
General solution: x = 179 (mod 1009)
>>> print_solutions(98, 105, 1001)
Solving the congruence: 98x = 105 (mod 1001)
Particular solution: x = 93
General solution: x = 93 (mod 143)
>>> print_solutions(98, 107, 1001)
Solving the congruence: 98x = 107 (mod 1001)
No solutions