解决大数的模块化线性同余
Solving modular linear congruences for large numbers
我正在寻找一种比我在 Whosebug 上找到的算法更好的算法来处理 4096 字节数字，我正在达到最大递归深度。
代码来自 stackoverlow post，我 copy/pasted 它但丢失了原来的 link:
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
这适用于较小的数字，例如：
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
但是我用更大的数字达到了最大递归深度。我提供的linear_congruence算法有没有更好的算法或者非递归算法？
根据 Eric Postpischil 的评论，我从维基百科条目中编写了伪代码，并利用此处的方法创建了一个非常快速的线性同余算法：http://gauss.math.luc.edu/greicius/Math201/Fall2012/Lectures/linear-congruences.article.pdf .
这对幂为 2-1 的 pow 很有效，可以得到答案。我正在研究如何抵消这个改变答案，并希望将它纳入这些答案也适用于这些答案，但现在，我有我需要的东西，因为我正在使用 pow 中 y 的 2 -1 的幂（ 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
它甚至可以即时处理 4096 字节的数字，这很棒：
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
感谢 Eric 指出这是什么，我使用 egcd 和上面 pdf 中的程序编写了一个非常快速的线性同余算法。如果任何 Whosebugers 需要一个快速算法，请将他们指向这个。我还了解到，当 pow(x,y,z) 的 y 偏离 2-1 的幂时，始终保持同余。我会进一步研究这个问题，看看是否存在偏移量变化以保持答案完整，如果找到，我会在未来跟进。
如果你有 Python 3.8 或更高版本，你可以用很少的代码行完成你需要的一切。
首先是一些数学：我假设你想求解 ax = b (mod m) 一个整数 x，给定整数 a、b 和 m.我还假设 m 是正数。
您首先需要计算的是 a 和 m 的最大公约数 g。有两种情况：
如果b不是g的倍数，则同余式无解（如果ax + my = b对于某些整数x和y，则 a 和 m 的公约数也必须是 b)
的公约数
如果b是g的倍数，那么同余正好等于(a/g)x = (b/g) (mod (m/g))。现在 a/g 和 m/g 互质，所以我们可以计算 a/g 模 m/g 的逆。将该逆乘以 b/g 得到一个解，并且可以通过将 m/g 的任意倍数添加到该解来获得通解。
Python 的 math module has had a gcd function since Python 3.5, and the built-in pow 函数可用于计算模逆，因为 Python 3.8.
综上所述，这是一些代码。首先是找到通用解决方案的函数，或者如果不存在解决方案则引发异常。如果成功，它 returns 两个整数。第一个给出了一个特定的解决方案；第二个给出提供通用解决方案的模数。
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
然后是一些驱动代码，演示如何使用上面的代码。
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})")
使用示例：
>>> 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
假设由于某种原因，您将 'attacking' 出现的线性同余方程 'empty'（无解）经常足以成为您算法的设计标准。
事实证明，您可以仅使用（有任何实际开销）剩余操作来回答该二进制问题 -
There exist solutions XOR There are no solutions
这可能对密码学有用；另见 abstract、
余数算术逻辑单元介绍
With Brief Computational Complexity Analysis
一旦确定存在解决方案，就可以使用反向替换
和 ALU 来确定解决方案。
此外，您将计算出 gcd(a,m) 并可以构造 Bézout 恒等式的系数
（如果你需要它们）。
以下是 python 结合了上述思想的程序；它计算存在的最小解并打印出 Bézout 的身份。
test_data = [ \
(32,12,82), \
(9,3,23), \
(17,41,73), \
(227,1,2011), \
(25,15,29), \
(2,22,71), \
(7,10,21), \
(124,58,900), \
(46, 12, 240), \
]
for lc in test_data:
LC = lc
back_sub_List = []
while True:
back_sub_List.append(LC)
n_mod_a = LC[2] % LC[0]
if n_mod_a == 0:
break
LC = (n_mod_a, -LC[1] % LC[0], LC[0])
gcd_of_a0_n0 = LC[0]
if LC[1] % LC[0] != 0:
print(f"No solution for {back_sub_List[0][0]}x = {back_sub_List[0][1]} (mod {back_sub_List[0][2]})")
else:
k = 0
for LC in back_sub_List[::-1]: # solve with back substitution
a,b,m = LC
k = (b + k*m) // a # optimize calculation since the remainder is zero?
print(f"The minimal solution for {back_sub_List[0][0]}x = {back_sub_List[0][1]} (mod {back_sub_List[0][2]}) is equal to {k}")
# get bezout
S = [1,0]
T = [0,1]
for LC in back_sub_List:
a,b,n = LC
q = n // a
s = S[0] - q * S[1]
S = [S[1], s]
t = T[0] - q * T[1]
T = [T[1], t]
print(f" Bézout's identity: ({S[0]})({lc[2]}) + ({T[0]})({lc[0]}) = {gcd_of_a0_n0}")
程序输出
The minimal solution for 32x = 12 (mod 82) is equal to 26
Bézout's identity: (-7)(82) + (18)(32) = 2
The minimal solution for 9x = 3 (mod 23) is equal to 8
Bézout's identity: (2)(23) + (-5)(9) = 1
The minimal solution for 17x = 41 (mod 73) is equal to 11
Bézout's identity: (7)(73) + (-30)(17) = 1
The minimal solution for 227x = 1 (mod 2011) is equal to 1320
Bézout's identity: (78)(2011) + (-691)(227) = 1
The minimal solution for 25x = 15 (mod 29) is equal to 18
Bézout's identity: (-6)(29) + (7)(25) = 1
The minimal solution for 2x = 22 (mod 71) is equal to 11
Bézout's identity: (1)(71) + (-35)(2) = 1
No solution for 7x = 10 (mod 21)
Bézout's identity: (0)(21) + (1)(7) = 7
No solution for 124x = 58 (mod 900)
Bézout's identity: (4)(900) + (-29)(124) = 4
The minimal solution for 46x = 12 (mod 240) is equal to 42
Bézout's identity: (-9)(240) + (47)(46) = 2
我正在寻找一种比我在 Whosebug 上找到的算法更好的算法来处理 4096 字节数字，我正在达到最大递归深度。
代码来自 stackoverlow post，我 copy/pasted 它但丢失了原来的 link:
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
这适用于较小的数字，例如：
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
但是我用更大的数字达到了最大递归深度。我提供的linear_congruence算法有没有更好的算法或者非递归算法？
根据 Eric Postpischil 的评论，我从维基百科条目中编写了伪代码，并利用此处的方法创建了一个非常快速的线性同余算法：http://gauss.math.luc.edu/greicius/Math201/Fall2012/Lectures/linear-congruences.article.pdf .
这对幂为 2-1 的 pow 很有效，可以得到答案。我正在研究如何抵消这个改变答案，并希望将它纳入这些答案也适用于这些答案，但现在，我有我需要的东西，因为我正在使用 pow 中 y 的 2 -1 的幂（ 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
它甚至可以即时处理 4096 字节的数字，这很棒：
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
感谢 Eric 指出这是什么，我使用 egcd 和上面 pdf 中的程序编写了一个非常快速的线性同余算法。如果任何 Whosebugers 需要一个快速算法，请将他们指向这个。我还了解到，当 pow(x,y,z) 的 y 偏离 2-1 的幂时，始终保持同余。我会进一步研究这个问题，看看是否存在偏移量变化以保持答案完整，如果找到，我会在未来跟进。
如果你有 Python 3.8 或更高版本，你可以用很少的代码行完成你需要的一切。
首先是一些数学：我假设你想求解 ax = b (mod m) 一个整数 x，给定整数 a、b 和 m.我还假设 m 是正数。
您首先需要计算的是 a 和 m 的最大公约数 g。有两种情况：
如果
的公约数b不是g的倍数，则同余式无解（如果ax + my = b对于某些整数x和y，则a和m的公约数也必须是b)如果
b是g的倍数，那么同余正好等于(a/g)x = (b/g) (mod (m/g))。现在a/g和m/g互质，所以我们可以计算a/g模m/g的逆。将该逆乘以b/g得到一个解，并且可以通过将m/g的任意倍数添加到该解来获得通解。
Python 的 math module has had a gcd function since Python 3.5, and the built-in pow 函数可用于计算模逆，因为 Python 3.8.
综上所述，这是一些代码。首先是找到通用解决方案的函数，或者如果不存在解决方案则引发异常。如果成功，它 returns 两个整数。第一个给出了一个特定的解决方案；第二个给出提供通用解决方案的模数。
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
然后是一些驱动代码，演示如何使用上面的代码。
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})")
使用示例：
>>> 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
假设由于某种原因，您将 'attacking' 出现的线性同余方程 'empty'（无解）经常足以成为您算法的设计标准。
事实证明，您可以仅使用（有任何实际开销）剩余操作来回答该二进制问题 -
There exist solutions XOR There are no solutions
这可能对密码学有用；另见 abstract、
余数算术逻辑单元介绍
With Brief Computational Complexity Analysis
一旦确定存在解决方案，就可以使用反向替换
和 ALU 来确定解决方案。
此外，您将计算出 gcd(a,m) 并可以构造 Bézout 恒等式的系数
（如果你需要它们）。
以下是 python 结合了上述思想的程序；它计算存在的最小解并打印出 Bézout 的身份。
test_data = [ \
(32,12,82), \
(9,3,23), \
(17,41,73), \
(227,1,2011), \
(25,15,29), \
(2,22,71), \
(7,10,21), \
(124,58,900), \
(46, 12, 240), \
]
for lc in test_data:
LC = lc
back_sub_List = []
while True:
back_sub_List.append(LC)
n_mod_a = LC[2] % LC[0]
if n_mod_a == 0:
break
LC = (n_mod_a, -LC[1] % LC[0], LC[0])
gcd_of_a0_n0 = LC[0]
if LC[1] % LC[0] != 0:
print(f"No solution for {back_sub_List[0][0]}x = {back_sub_List[0][1]} (mod {back_sub_List[0][2]})")
else:
k = 0
for LC in back_sub_List[::-1]: # solve with back substitution
a,b,m = LC
k = (b + k*m) // a # optimize calculation since the remainder is zero?
print(f"The minimal solution for {back_sub_List[0][0]}x = {back_sub_List[0][1]} (mod {back_sub_List[0][2]}) is equal to {k}")
# get bezout
S = [1,0]
T = [0,1]
for LC in back_sub_List:
a,b,n = LC
q = n // a
s = S[0] - q * S[1]
S = [S[1], s]
t = T[0] - q * T[1]
T = [T[1], t]
print(f" Bézout's identity: ({S[0]})({lc[2]}) + ({T[0]})({lc[0]}) = {gcd_of_a0_n0}")
程序输出
The minimal solution for 32x = 12 (mod 82) is equal to 26
Bézout's identity: (-7)(82) + (18)(32) = 2
The minimal solution for 9x = 3 (mod 23) is equal to 8
Bézout's identity: (2)(23) + (-5)(9) = 1
The minimal solution for 17x = 41 (mod 73) is equal to 11
Bézout's identity: (7)(73) + (-30)(17) = 1
The minimal solution for 227x = 1 (mod 2011) is equal to 1320
Bézout's identity: (78)(2011) + (-691)(227) = 1
The minimal solution for 25x = 15 (mod 29) is equal to 18
Bézout's identity: (-6)(29) + (7)(25) = 1
The minimal solution for 2x = 22 (mod 71) is equal to 11
Bézout's identity: (1)(71) + (-35)(2) = 1
No solution for 7x = 10 (mod 21)
Bézout's identity: (0)(21) + (1)(7) = 7
No solution for 124x = 58 (mod 900)
Bézout's identity: (4)(900) + (-29)(124) = 4
The minimal solution for 46x = 12 (mod 240) is equal to 42
Bézout's identity: (-9)(240) + (47)(46) = 2