任意 X 和 Y 的最快整数搜索

Fastest Integer Search for Any Arbitrary X and Y

随着离散数学的应用,python 解决这个问题最快的算法是什么:

用等式ax + by = d,其中a, bd 是用户输入,搜索 [=39= 的整数值对]x 和 y 在 |L, R| 范围内(包括L和R)满足方程

L 和 R 也是用户输入。

按升序输出 x 和 y 的所有可能值。如果没有可能的对,打印 none

案例 1:

a = 1
b = 5
d = 40
L = 0
R = 10

Result:
0 8
5 7
10 6

案例二:

a = 14
b = 91
d = 53
L = 1
R = 100

Result:
none

这是我的代码,但我相信有一种 更快的搜索方法。这样效率太低了。

a = int(input())
b = int(input())
d = int(input())
L = int(input())
R = int(input())

isNone = True
for x in range(L, R+1):
    for y in range(L, R+1):
        if (a*x) + (b*y) == d:
            print(x, y)
            isNone = False
if isNone: print("none")

可以有一个算法O(1)吗?最快的方法是什么?

我想你想要申请的身份如下(直接来自维基百科,虽然在大多数离散数学文本中应该可以找到类似的措辞,或者可以自己证明):

The simplest linear Diophantine equation takes the form ax + by = c, where a, b and c are given integers. The solutions are described by the following theorem:

This Diophantine equation has a solution (where x and y are integers) if and only if c is a multiple of the greatest common divisor of a and b. Moreover, if (x, y) is a solution, then the other solutions have the form (x + kv, y − ku), where k is an arbitrary integer, and u and v are the quotients of a and b (respectively) by the greatest common divisor of a and b.

这几乎立即回答了这个问题,特别是因为它的标准证明使用了一种叫做 euclidean 算法 的东西。为了简单起见,我们将执行以下操作:

  1. 正向使用欧氏算法求g = gcd(a, b).
  2. 通过欧几里德算法反向求解找到 _x, _y 使得 _x*a + _y*b == g.
  3. 如果d不是g的倍数那么不可能有任何解决方案,所以早点退出。
  4. 否则,x, y = _x*(d//g), _y*(d//g)是一个可能的解决方案。用它来找到所需范围内的所有解决方案。
def gcd(a, b):
    # forward euclidean algorithm
    q, r, x, qs = None, b, a, []
    while x%r:
        (q, r), x = divmod(x, r), r
        qs.append(q)

    # save the gcd for later
    g = r

    # backsolve euclidean algorithm
    if not qs:
        return 1, 1-a//b, g
    theta, omega = 1, -qs[-1]
    for q in reversed(qs[:-1]):
        theta, omega = omega, theta - q * omega
    
    # theta * a + omega * b == g
    # g might be negative, but we don't care about a canonical solution
    return theta, omega, g

def idivide_zero(a, b):
    # integer division a/b, round toward 0 instead of round down
    q = a // b
    if q < 0 and b*q != a:
        q += 1
    return q

def bounded_solutions(a, b, d, L, R):
    _x, _y, g = gcd(a, b)
    if d%g:
        return

    # a*x + b*y == d
    x, y = _x*(d//g), _y*(d//g)

    # solutions are of the form (x+k*v, y-k*u)
    u, v = a//g, b//g

    # The next trick is to find all solutions in [L, R].
    # Basically, we need L <= x+k*v <= R and L <= y-k*u <= R.
    # Note that valid choices of k exist in a contiguous interval, so
    # we only have to find the lower and upper bounds to be able to
    # quickly enumerate all options.
    xb = sorted(idivide_zero(b-x, v) for b in (L, R))
    yb = sorted(idivide_zero(y-b, u) for b in (L, R))
    m, M = min(xb[0], yb[0]), max(xb[1], yb[1])
    for k in range(m, M+1):
        yield x+k*v, y-k*u

a = int(input())
b = int(input())
d = int(input())
L = int(input())
R = int(input())

empty = True
for x, y in bounded_solutions(a, b, d, L, R):
    print(x, y)
    empty = False
if empty:
    print('none')

代码未经测试。精神上是正确的,但可能还有一些调试遗留。