c++ - 方程 `a + bx = c + dy` 的积分解

Question

在等式 a + bx = c + dy 中，所有变量都是整数。 a、b、c 和 d 是已知的。如何找到 x 和 y 的积​​分解？如果我的想法是正确的，将有无限多的解决方案，由 b 和 d 的最小公倍数分隔，但我只需要一个解决方案，并且我可以计算其余部分。这是一个例子:

a = 2
b = 3
c = 4
d = 5
a + bx: (2, 5,  8, 11, 14)
c + dy: (4, 9, 14, 19, 24)

a + bx intersects c + dy at 14, so:
x = 4
y = 2

现在，我正在遍历 x 的整数值，直到找到 y 的整数值(伪代码):

function integral_solution(int a, int b, int c, int d) {
    // a + bx == c + dy
    // (a + bx - c) / d == y

    // Some parameters may have no integral solution,
    // for example if b == d and (a - c) % b != 0
    // If we don't have a solution before x == d, there is none.

    int x = 0;
    while (x < d)
    {
        if ((a + bx - c) % d == 0)
        {
            return [x, (a + bx - c) / d];
        }
        x++;
    }
    return false;
}

我觉得有更好的方法来做到这一点。有什么办法不用循环就可以找到 x 和 y 吗？我正在使用 C++，如果这很重要的话。

Answer 1

Linear Diophantine方程采用 ax + by = c 的形式。如果 c 是 a 和 b 的最大公约数，这意味着 a=z'c 和 b=z''c 那么这是 Bézout's identity形式的

a=z' 和 b=z'' 方程有无限多解。因此，除了尝试搜索方法之外，您还可以检查c 是否是a 和b 的最大公约数 (GCD) (在您的情况下，这转化为 bx - dy = c - a)

如果 a 和 b 确实是 c 的倍数，那么 x 和 y 可以使用 extended Euclidean algorithm 计算它找到满足 Bézout 身份的整数 x 和 y(其中一个通常为负数)

你的答案是:

a = k*x,b = k*y,c - a = k * gcd(a,b) 对于任何整数 k。

(作为旁注:这也适用于任何其他 Euclidean domain ，即多项式环和每个欧几里德域都是 unique factorization domain )。您可以使用迭代法找到这些解决方案:

---

迭代法

通过对同类项进行展开分组的常规代数运算(引用前文wikipedia article的最后一节)，迭代法得到如下算法:

• 1。应用欧氏算法，令qn(n从1开始)为除法中商的有限列表。
• 2。分别将x0、x1初始化为1、0，将y0、y1初始化为0,1。
  • 2.1 然后对于每个 i 只要定义了 qi，
  • 2.2 计算xi+1 = xi−1 − qixi
  • 2.3 计算yi+1 = yi−1 − qiyi
  • 2.4 i 加 1 后重复上述操作。
• 3。答案是 xn 和 yn 的倒数第二个。

伪代码:

function extended_gcd(a, b)
    x := 0    lastx := 1
    y := 1    lasty := 0
    while b ≠ 0
        quotient := a div b
        (a, b) := (b, a mod b)
        (x, lastx) := (lastx - quotient*x, x)
        (y, lasty) := (lasty - quotient*y, y)       
    return (lastx, lasty)

因此我编写了示例算法，该算法使用欧几里德算法 迭代方法计算最大公约数 非负a 和b (对于负值 - these 需要额外的步骤)，它返回 GCD 并将 x 和 y 的解决方案存储在通过引用传递给它的变量中:

int gcd_iterative(int a, int b, int& x, int& y) {
    int c;
    std::vector<int> r, q, x_coeff, y_coeff;
    x_coeff.push_back(1); y_coeff.push_back(0);
    x_coeff.push_back(0); y_coeff.push_back(1);

    if ( b == 0 ) return a;
    while ( b != 0 ) {
            c = b;
            q.push_back(a/b);
            r.push_back(b = a % b);
            a = c;
            x_coeff.push_back( *(x_coeff.end()-2) -(q.back())*x_coeff.back());
            y_coeff.push_back( *(y_coeff.end()-2) -(q.back())*y_coeff.back());
    }
    if(r.size()==1) {
        x = x_coeff.back();
        y = y_coeff.back();
    } else {
        x = *(x_coeff.end()-2);
        y = *(y_coeff.end()-2);
    }
    std::vector<int>::iterator it;
    std::cout << "r: ";
    for(it = r.begin(); it != r.end(); it++) { std::cout << *it << "," ; }
    std::cout << "\nq: ";
    for(it = q.begin(); it != q.end(); it++) { std::cout << *it << "," ; }
    std::cout << "\nx: ";
    for(it = x_coeff.begin(); it != x_coeff.end(); it++){ std::cout << *it<<",";}
    std::cout << "\ny: ";
    for(it = y_coeff.begin(); it != y_coeff.end(); it++){ std::cout << *it<<",";}
    return a;
} 

通过传递给它一个 example从维基百科 a = 120 和 b = 23 我们得到:

int main(int argc, char** argv) {   
    // 120x + 23y = gcd(120,23)
    int x_solution, y_solution;
    int greatestCommonDivisor =  gcd_iterative(120, 23, x_solution, y_solution);
    return 0;
}

r: 5,3,2,1,0,

q: 5,4,1,1,2,

x: 1,0,1,-4,5,-9,23,

y: 0,1,-5,21,-26,47,-120,

根据这个例子给定的表格是什么: