自然数域ax+by=c求解：扩展欧几里得算法适配性问询

我是一名离散数学方向的业余爱好者，本问题不属于课后作业，仅为个人居家研究时遇到的疑问。
该问题形式类似贝祖恒等式（Bezout's identity），但求解域限定为自然数，等式右端为任意给定常数c。在给定自然数a、b、c的前提下，需要找出所有满足等式的自然数解对(x, y)。
整数域下的贝祖恒等式存在无穷多组整数解，二者的结构相似性意味着扩展欧几里得算法（extended Euclidian algorithm, EEA）可用于本问题求解。以下是从公开资料改编、目前可正常运行的两个EEA实现（递归版与迭代版）：

typedef long int Int;

#ifdef RECURSIVE_EEA
Int // returns the GCD of a and b and finds x and y 
// such that ax + by == GCD(a,b), recursively
eea(Int a, Int b, Int &x, Int &y) {
    if (0==a) {
        x = 0;
        y = 1;
        return b;
    }
    Int x1; x1=0;
    Int y1; y1=0;
    Int gcd = eea(b%a, a, x1, y1);
    x = y1 - b/a*x1;
    y = x1;
    return gcd;
}
#endif

#ifdef ITERATIVE_EEA
Int // returns the GCD of a and b and finds x and y 
// such that ax + by == GCD(a,b), iteratively
eea(Int a, Int b, Int &x, Int &y) {
    x = 0;
    y = 1; 
    Int u; u=1;
    Int v; v=0;         // does this need initialising?
    Int q;              // quotient
    Int r;              // remainder
    Int m;
    Int n;
    while (0!=a) {      
        q = b/a;        // quotient
        r = b%a;        // remainder
        m = x - u*q;    // ??  what are the invariants?
        n = y - v*q;    // ??  When does this overflow?
        b = a;          // A candidate for the gcd - a's last nonzero value.
        a = r;          // a becomes the remainder - it shrinks each time.
                        // When a hits zero, the u and v that are written out
                        // are final values and the gcd is a's previous value.
        x = u;          // Here we have u and v shuffling values out
        y = v;          // via x and y. If a has gone to zero, they're final.
        u = m;          // ... and getting new values
        v = n;          // from m and n
    }
    return b;
}
#endif

核心疑问：这两段代码能否改造适配上述自然数域方程的求解需求？是否存在其他更优的求解思路？

解答 #

现有EEA代码的改造方案 #

你写的递归、迭代两个版本EEA都可以直接改造适配需求，不需要修改核心EEA逻辑，只需要在外层增加解的缩放、范围筛选步骤即可：

• 第一步：可解性校验
  调用EEA计算g = eea(a,b,x0,y0)，得到a和b的最大公约数g，以及满足a*x0 + b*y0 = g的整数特解(x0,y0)。如果c不能被g整除，方程不存在任何自然数解，直接返回空结果。
• 第二步：缩放得到目标方程特解
  令缩放系数k = c / g，得到整数域下ax + by = c的一个特解：

  x_s = x0 * k
  y_s = y0 * k
  

• 第三步：推导整数通解，筛选自然数范围
  该方程的所有整数解都可以表示为：

  x = x_s + t * (b/g)
  y = y_s - t * (a/g)
  

  其中t为任意整数。加上自然数约束（若定义自然数为非负整数则取x≥0、y≥0；若定义为正整数则取x≥1、y≥1），解出t的合法取值区间：
  • 由x的约束得到t的下界：t_min = ceil( (约束下界 - x_s) * g / b )
  • 由y的约束得到t的上界：t_max = floor( (y_s - 约束下界) * g / a )
    如果t_min > t_max，则无自然数解；否则所有t取t_min到t_max之间的整数时，对应的(x,y)就是全部自然数解。

针对你代码里的注释疑问：迭代版中u和v必须初始化，你当前写的u=1; v=0;是完全正确的；循环的不变量是每轮迭代都满足a*u + b*v == r、a*x + b*y == b（即当前迭代的两个余数都可以被a、b线性表示）；溢出仅在a、b取值接近long int上限时才会出现，常规研究场景下不需要额外处理。两段代码的核心逻辑都没有问题，计算出的gcd和特解都是准确的。

更优的轻量求解思路 #

针对自然数域的这个特定问题，有比通用EEA实现更简单、更不容易出错的方案，时间复杂度和EEA一致：

1. 先做gcd校验：如果c不能被g = gcd(a,b)整除，直接返回无解；否则将方程两边同除以g，得到等价简化方程a'x + b'y = c'，其中a'=a/g、b'=b/g、c'=c/g，此时a'和b'互质。
2. 求解模运算下的最小x：因为a'和b'互质，a'在模b'下存在乘法逆元，可直接算出满足a' * x0 ≡ c' mod b'的最小非负整数x0，对应的y0 = (c' - a'*x0)/b'，这就是x最小的一组特解。
3. 生成所有解：所有自然数解的x取值为x = x0 + k*b'（k为非负整数），直到x超过c'/a'为止，对应的y值为y = (c' - a'*x)/b'，全部为自然数。
   如果不想实现模逆元逻辑，也可以直接遍历x从0到c'/a'，判断(c' - a'*x)是否能被b'整除，在a、b数值不大的场景下代码量极小，运行效率也足够。

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