求平面Ax+By+Cz-D=0上非负整数点的高效算法咨询

问题描述 #

给定平面方程 Ax + By + Cz = D（满足D远大于A、B、C，且GCD(A,B,C)=1），需要找出所有满足x、y、z为非负整数的点(x,y,z)。原暴力解法通过遍历x、y的可能范围计算z并验证，当D增大时效率极低，需更高效的算法或优化方案。

原暴力Python代码（注：原代码存在变量名错误，已修正明显问题）：

points = []
a, b, c, d =  # 需替换为实际参数
for x in range(d // a + 1):
    remaining = d - a * x
    if remaining < 0:
        continue
    for y in range(remaining // b + 1):
        z_remain = remaining - b * y
        if z_remain % c == 0:
            z = z_remain // c
            points.append((x, y, z))

原代码存在变量名错误：循环中用了z作为x的范围除数，应该是a；内层循环用y作为除数，应该是b；最后append的是(x,y,x)，应该是(x,y,z)，以上已修正。

优化思路与高效算法 #

1. 降维优化：固定一维，转化为二元一次不定方程求解 #

因为GCD(A,B,C)=1，我们可以先固定其中一个变量（比如x），将方程转化为 By + Cz = D - Ax，此时只需确保D - Ax ≥ 0，然后求解这个二元一次不定方程的非负整数解(y,z)。

对于二元一次不定方程 Py + Qz = R（这里P=B, Q=C, R=D-Ax），求解步骤：

• 首先判断R是否≥0，且GCD(P,Q)能整除R（因为原GCD(A,B,C)=1，当R=D-Ax时，GCD(B,C)必然整除R，否则该x对应的方程无解）。
• 用扩展欧几里得算法找到一组特解(y₀,z₀)，然后根据通解公式生成所有非负解：
  通解为 y = y₀ + (Q/g) * t，z = z₀ - (P/g) * t，其中g=GCD(P,Q)，t为整数。
• 找出所有t的取值范围，使得y≥0且z≥0，即可得到该x对应的所有(y,z)解。

2. 减少遍历维度 #

相比暴力法的两层循环，这种方法只需要遍历单个变量的可能值，每个变量对应的求解是O(1)级别的（只需计算t的范围），整体时间复杂度从O(D²/(A*B))降低到O(D/max(A,B,C))，效率提升明显。

3. 进一步优化：选择遍历范围最小的维度 #

不要固定x，而是选择A、B、C中最大的那个对应的变量来遍历，比如如果A是最大的，那么x的范围是0到D//A，这个范围最小，遍历次数最少。

实现示例（Python） #

import math

def extended_gcd(p, q):
    # 扩展欧几里得算法，返回(g, x, y)，满足p*x + q*y = g
    if q == 0:
        return (p, 1, 0)
    else:
        g, x, y = extended_gcd(q, p % q)
        return (g, y, x - (p // q) * y)

def solve_equation(a, b, c, d):
    points = []
    # 选择遍历系数最大的变量，减少循环次数
    coeffs = [(a, 'x'), (b, 'y'), (c, 'z')]
    coeffs.sort(reverse=True, key=lambda x: x[0])
    max_coeff, var = coeffs[0]
    
    if var == 'x':
        for x in range(0, d // a + 1):
            r = d - a * x
            if r < 0:
                continue
            g, y0, z0 = extended_gcd(b, c)
            if r % g != 0:
                continue
            # 调整特解到满足方程的形式
            k = r // g
            y0 *= k
            z0 *= k
            # 通解参数t的范围
            step_y = c // g
            step_z = b // g
            # 计算t的下限：y ≥ 0
            t_min = math.ceil(-y0 / step_y) if step_y != 0 else 0
            # 计算t的上限：z ≥ 0
            t_max = math.floor(z0 / step_z) if step_z != 0 else 0
            # 遍历所有合法的t
            for t in range(t_min, t_max + 1):
                y = y0 + step_y * t
                z = z0 - step_z * t
                if y >= 0 and z >= 0:
                    points.append((x, y, z))
    elif var == 'y':
        # 遍历y的逻辑，类似x
        for y in range(0, d // b + 1):
            r = d - b * y
            if r < 0:
                continue
            g, x0, z0 = extended_gcd(a, c)
            if r % g != 0:
                continue
            k = r // g
            x0 *= k
            z0 *= k
            step_x = c // g
            step_z = a // g
            t_min = math.ceil(-x0 / step_x) if step_x != 0 else 0
            t_max = math.floor(z0 / step_z) if step_z != 0 else 0
            for t in range(t_min, t_max + 1):
                x = x0 + step_x * t
                z = z0 - step_z * t
                if x >= 0 and z >= 0:
                    points.append((x, y, z))
    else:
        # 遍历z的逻辑，类似x
        for z in range(0, d // c + 1):
            r = d - c * z
            if r < 0:
                continue
            g, x0, y0 = extended_gcd(a, b)
            if r % g != 0:
                continue
            k = r // g
            x0 *= k
            y0 *= k
            step_x = b // g
            step_y = a // g
            t_min = math.ceil(-x0 / step_x) if step_x != 0 else 0
            t_max = math.floor(y0 / step_y) if step_y != 0 else 0
            for t in range(t_min, t_max + 1):
                x = x0 + step_x * t
                y = y0 - step_y * t
                if x >= 0 and y >= 0:
                    points.append((x, y, z))
    return points

关键说明 #

• 扩展欧几里得算法是求解不定方程的核心，用于找到方程的特解。
• 选择遍历系数最大的变量，能最大程度减少循环次数，尤其是当D远大于系数时，效果更明显。
• 通解的t范围计算需要注意整数边界，使用math.ceil和math.floor确保覆盖所有合法的整数解。

内容的提问来源于stack exchange，提问作者Ilikecoding