解决大量的模块化线性同余 [英] Solving modular linear congruences for large numbers

问题描述

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

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.

stackoverlow帖子中的代码，我复制/粘贴了它，但丢失了原始链接:

Code from stackoverlow post, i copy/pasted it but lost the original 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

这对于较小的数字很好用，例如:

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

但是我用更大的数字达到了最大递归深度.我提供的linear_congruence算法是否有更好的算法或非递归算法?

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?

基于Eric Postpischil的评论，我从Wikipedia条目中编写了伪代码，并使用以下方法创建了一种非常快速的线性一致算法:

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 .

这对于功率为2-1的战俘非常有效，以获得答案.我正在研究如何抵消这个改变答案，并希望将它纳入这些答案，但现在，我有我需要的东西，因为我在 pow(x，y，z):

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

它甚至可以立即处理4096个字节的数字，这很不错:

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

感谢Eric指出这是什么，我写了一个非常快速的线性一致算法，它使用egcd和上面pdf中的过程.如果有任何stackoverflowers需要快速算法，请将其指向该算法.我还了解到，当pow(x，y，z)的y取2-1的幂时，始终保持一致.我将对此进行进一步调查，以查看是否存在偏移量更改以保持答案的完整性，并且如果发现该错误，将在以后进行跟进.

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.

推荐答案

如果您使用的是Python 3.8或更高版本，则只需很少的几行代码即可完成所需的一切.

If you have Python 3.8 or later, you can do everything you need to with a very small number of lines of code.

首先是一些数学运算:我假设您要为整数 x 求解 ax = b(mod m)，对于给定的整数 a ， b 和 m .我还假设 m 为正.

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.

您需要计算的第一件事是 a 和 m 的最大公约数 g .有两种情况:

The first thing you need to compute is the greatest common divisor g of a and m. There are two cases:

• 如果 b 不是 g 的倍数，则全等式无解(如果 ax + my = b 用于一些整数 x 和 y ，那么 a 和 m 的任何常见除数也必须是的除数> b )

• 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)

如果 b 是 g 的倍数，则等价性完全等同于(a/g)x =(b/g)(mod(m/g)).现在 a/g 和 m/g 是相对质数，因此我们可以对 a/g 取模 m/g.将该逆乘以 b/g 即可得出一个解决方案，可以通过将任意 m/g 倍数添加到该解决方案中来获得一般解决方案.

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的 数学 模块具有 gcd 自Python 3.5以来的功能，以及内置的 pow 函数可用于自Python 3.8起计算模块化逆.

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})")

示例用法:

>>> 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

这篇关于解决大量的模块化线性同余的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！