求解大数的模线性同余

我正在寻找一种比我在 stackoverflow 上找到的算法更好的算法来处理 4096 字节数，我正在达到最大递归深度。

来自 stackoverflow 帖子的代码，我复制/粘贴了它，但丢失了原始链接：

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

但我用更大的数字达到了最大递归深度。我提供的线性同余算法是否有更好的算法或非递归算法？

根据 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 中的过程编写了一个非常快速的线性同余算法。如果任何 stackoverflowers 需要快速算法，请向他们指出这个。我还了解到，当 pow(x,y,z) 的 y 具有 2-1 的幂时，始终保持同余。我将进一步研究这一点，看看是否存在偏移变化以保持答案完整，如果发现，我将在将来跟进。

如果您有 Python 3.8 或更高版本，您可以用很少的代码行完成您需要的一切。

首先是一些数学：我假设你想解决ax = b (mod m)对于一个整数x，给定整数a, b and m。我也假设m是积极的。

您需要计算的第一件事是最大公约数g of a and m。有两种情况：

• if b不是的倍数g，则同余无解（如果ax + my = b对于某些整数x and y，那么任何公约数a and m也必须是除数b)

• if b is的倍数g，那么同余式完全等价于(a/g)x = (b/g) (mod (m/g)). Now a/g and m/g互质，因此我们可以计算其倒数a/g modulo m/g。将该倒数乘以b/g给出一个解，将任意倍数相加即可得到通解m/g到该解决方案。

蟒蛇的math模块有一个gcd自Python 3.5以来的函数，以及内置的pow自 Python 3.8 起，函数可用于计算模逆。

将它们放在一起，这是一些代码。首先是一个找到通用解决方案的函数，或者如果不存在解决方案则引发异常。如果成功，则返回两个整数。第一个给出了特定的解决方案；第二个给出提供通解的模数。

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