基于离散数学：Python求解ax+by=d整数解的最快算法问询

阿华AIGC实验室

2026-5-9

线性丢番图方程区间解的高效求解方法

首先直接回答你的问题：不存在严格意义上的O(1)算法——因为如果区间内有k个解，输出这些解至少需要O(k)的时间。但我们可以实现**时间复杂度为O(log(min(|a|,|b|)) + k)**的算法，这比你当前的O((R-L+1)^2)暴力解法快几个数量级，尤其是当R-L很大时（比如1e5甚至更大）。

核心思路：基于数论的精确求解

线性丢番图方程ax + by = d的解的存在性和通解都有明确的数论结论，我们可以利用这些结论直接计算出区间[L,R]内的所有解，而不需要暴力枚举：

1. 解的存在性判断

方程有整数解的充要条件是：gcd(a, b) 整除 d。记g = gcd(a, b)，如果d % g != 0，直接输出none即可。

2. 简化方程

如果方程有解，我们可以将方程两边同时除以g，得到简化后的方程：a'x + b'y = d'，其中：

a' = a // g
b' = b // g
d' = d // g
此时gcd(a', b') = 1，这一步是为了方便后续求特解和通解。

3. 求特解

我们需要找到简化方程a'x + b'y = d'的一组特解(x0, y0)。这可以通过扩展欧几里得算法实现：先找到满足a'x + b'y = 1的解，再将解乘以d'得到目标特解。

4. 通解推导

因为gcd(a', b') = 1，简化方程的通解可以表示为：

x = x0 + t * b'
y = y0 - t * a'

其中t是任意整数。（注：通解的形式可能因特解的计算方式略有不同，核心是保证代入原方程成立，调整t的系数符号即可）

5. 确定t的有效范围

我们需要找到所有整数t，使得x ∈ [L, R]且y ∈ [L, R]：

对于x的约束：L ≤ x0 + t*b' ≤ R
对于y的约束：L ≤ y0 - t*a' ≤ R
解这两个不等式，得到t的取值范围[t_low, t_high]（注意当系数为负数时，不等式方向要反转）。如果t_low > t_high，说明区间内没有解，输出none。

6. 生成并输出解

遍历t从t_low到t_high的所有整数，计算对应的x和y，按x升序输出（因为x随t的变化方向由b'的符号决定，调整遍历顺序即可保证升序）。

Python实现代码

import math

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

def solve_linear_diophantine(a, b, d, L, R):
    g, x0, y0 = extended_gcd(a, b)
    # 判断解是否存在
    if d % g != 0:
        print("none")
        return
    
    # 简化方程
    a_prime = a // g
    b_prime = b // g
    d_prime = d // g
    
    # 计算特解
    x0 *= d_prime
    y0 *= d_prime
    
    # 确定t的范围：x = x0 + t*b_prime ∈ [L, R]
    def get_t_range(val, coeff, lower, upper):
        # 解不等式 lower ≤ val + t*coeff ≤ upper，返回t的整数范围[low, high]，无解则返回(1,0)
        if coeff == 0:
            # 系数为0，检查val是否在区间内
            if lower <= val <= upper:
                # t可以是任意整数，但需要结合另一个约束
                return (-float('inf'), float('inf'))
            else:
                return (1, 0)
        elif coeff > 0:
            t_low = math.ceil((lower - val) / coeff)
            t_high = math.floor((upper - val) / coeff)
        else:
            # 系数为负，不等式方向反转
            t_low = math.floor((upper - val) / coeff)
            t_high = math.ceil((lower - val) / coeff)
        return (t_low, t_high)
    
    t_x_low, t_x_high = get_t_range(x0, b_prime, L, R)
    t_y_low, t_y_high = get_t_range(y0, -a_prime, L, R)
    
    # 取交集
    t_low = max(t_x_low, t_y_low)
    t_high = min(t_x_high, t_y_high)
    
    if t_low > t_high:
        print("none")
        return
    
    # 按x升序输出解：如果b_prime是正的，t从低到高x递增；如果b_prime是负的，t从高到低x递增
    solutions = []
    if b_prime >= 0:
        t_start, t_end, step = t_low, t_high, 1
    else:
        t_start, t_end, step = t_high, t_low, -1
    
    for t in range(t_start, t_end + step, step):
        x = x0 + t * b_prime
        y = y0 - t * a_prime
        solutions.append((x, y))
    
    # 输出所有解，按要求格式
    if not solutions:
        print("none")
    else:
        print(' '.join(f"{x} {y}" for x, y in solutions))

# 读取输入
a = int(input())
b = int(input())
d = int(input())
L = int(input())
R = int(input())

solve_linear_diophantine(a, b, d, L, R)