数学/数论专题-学习笔记:扩展欧几里得(exgcd)

1. 前言

扩展欧几里得(exgcd),是在欧几里得算法基础上求解任意形如 \(ax+by=c\) 的二元一次方程的一组特解的一种算法。

在往下看之前,您只需要知道如何使用欧几里得算法求 \(\gcd(a,b)\)

不知道也没关系,式子在这里:

\[\gcd(a,b)=\gcd(b, a \bmod b) \]

证明网络上面有很多,可以自行查找。

那么 \(\text{exgcd}\) 就是其扩展算法。

2. 详解

2.1 exgcd 过程

例题 1:给出 \(a,b \in N_+,d=\gcd(a,b)\),求方程 \(ax+by=d\) 的一组整数解。

这就是 \(\text{exgcd}\) 要解决的问题。

接下来给出证明 \(ax+by=d\) 一定有整数解的证明过程,该过程同时也给出了 \(\text{exgcd}\) 的求解过程。

考虑普通 \(\gcd\) 的求解过程:\(\gcd(a,b)=\gcd(b,a \bmod b)\)

假设我们现在已经知道了一组解 \(x',y'\) 是方程 \(bx + (a \bmod b)y=d\) 的一组特解,那么如何求出 \(ax+by=d\) 的解呢?

首先我们知道这个式子成立:$$a \bmod b=a - \left\lfloor\dfrac{a}{b}\right\rfloor \times b$$

那么把上面这个式子带入到 \(bx+(a \bmod b)y=d\) 中,就有:

\[bx'+(a-\left\lfloor\dfrac{a}{b}\right\rfloor \times b)y'=d \]

拆掉括号:

\[bx' + ay' -\left\lfloor\dfrac{a}{b}\right\rfloor \times b \times y' = d \]

左边转化一下:

\[ay'+b(x'-\left\lfloor\dfrac{a}{b}\right\rfloor y')=d \]

上述方程跟下述方程是等价的:

\[ax+by=d \]

两者对比,得到:

\[\begin{cases}x=y'\\y=x'-\left\lfloor\dfrac{a}{b}\right\rfloor y'\end{cases} \]

于是我们成功的通过 \(bx'+(a \bmod b)y'=d\) 推出了 \(ax+by=d\) 的解。

那么为什么一定会有解呢?

考虑欧几里得算法求 \(\gcd\) 的步骤:\(\gcd(a,b)=\gcd(b,a \bmod b)\)

仿照上述过程求 \(\text{exgcd}\),那么初始状态是 \(a'x+b'y=d\),其中 \(a=d,b=0\)

因此 \(a=1,b=0\)

然后根据上述分析不断还原,最后一定能够还原出 \(ax+by=d\) 的解。

证毕。

上述过程给出了 \(ax+by=d\) 一定有解的证明,同时阐述了为什么有解。

同时根据上述过程应该明白为什么右边是 \(\gcd(a,b)=d\) 了吧~这跟 \(\gcd\) 的求法有很大关系。

Code:

void exgcd(int a, int b, LL &x, LL &y)//注意 x 跟 y 是引用
{
    if (b == 0) {x = 1; y = 0; return ;}
    exgcd(b, a % b, x, y);
    LL p = x; x = y; y = p - ((LL)a / b) * y;//特别注意存一下 x 的初始值
}

int main()
{
	a = Read(), b = Read(); LL x, y;
	exgcd(a, b, x, y); printf("%lld %lld\n", x, y);
}

2.2 一般情况

例题 2:P5656 【模板】二元一次不定方程 (exgcd)

这道题有一点恶心,主要是因为 0 不是正整数(

本题要求解决 \(ax+by=c\) 的正整数解的若干问题:判定无整数解,有无正整数解,有正整数解的情况下正整数解中 \(x,y\) 的最小最大值。

首先根据裴蜀定理,得知如果 \(\gcd(a,b)\not\mid c\) 则原方程无整数解。

否则设我们用 exgcd 求出方程 \(ax+by=\gcd(a,b)\) 的一组特解为 \(x',y'\),那么原方程的一组特解为 \(x_0=x'\times\dfrac{c}{\gcd(a,b)},y_0=y'\times\dfrac{c}{\gcd(a,b)}\),其实就是两边同乘 \(\dfrac{c}{\gcd(a,b)}\) 使得等式右边为 \(c\)

根据数学知识,在实数域上 \(x,y\) 的通解为:

\[\begin{cases}x=x_0+mb\\y=y_0-ma\end{cases},m\in R \]

但是要求 \(x,y\) 是整数,而为了覆盖所有符合要求的整数解则要求 \(mb,ma\) 尽量小,据此知道 \(m=k\times\dfrac{1}{\gcd(a,b)},k\in Z\) 可以使 \(ma,mb\) 尽量小,因此在整数域上 \(x,y\) 的通解为:

\[\begin{cases}x=x_0+k\times\dfrac{b}{\gcd(a,b)}\\y=y_0-k\times\dfrac{a}{\gcd(a,b)}\end{cases},k\in Z \]

接下来考虑按照例题要求判断有无正整数解,显然如果 \(x\) 越小那么 \(y\) 越大,因此我们可以考虑求出 \(x\) 的最小正整数值,注意到通解形式就是不断往 \(x_0\) 上加上或减去 \(\dfrac{b}{\gcd(a,b)}\),因此最小正整数值显然在 \([1,\dfrac{b}{\gcd(a,b)}]\) 内,因此我们直接将 \(x_0\)\(\dfrac{b}{\gcd(a,b)}\) 取模即可(特别注意 \(x=0\) 的情况要化成 \(\dfrac{b}{\gcd(a,b)}\),后面的所有取模同理,不再赘述),然后用 \(ax+by=c\) 算出 \(y\),此时在正整数解内有 \(x\) 最小 \(y\) 最大,但若此时 \(y\) 依然不是正整数那么就是无正整数解的情况,仿照求 \(x\) 最小值求出 \(y\) 最小值即可。

剩下的情况就是有正整数解的情况了,在之前我们已经求出了 \(x\) 最小 \(y\) 最大的解,类似的求出 \(x\) 最大 \(y\) 最小的解就得到了 \(x,y\) 最小值最大值,至于解的个数,由于相邻两个正整数解 \(x_i,x_{i+1}\) 差值就是 \(\dfrac{b}{\gcd(a,b)}\),因此解的个数就是 \(\dfrac{\max(x)-\min(x)}{\frac{b}{\gcd(a,b)}}+1\)(不理解的考虑等差数列求项数)。

写代码的时候如果要写 \(x\)\(p\) 取模可以这么写:x = (x % p + p) % p,然后注意一下全部开 long long 以及 \(x,y\) 最小值不能有 0 即可。

Code:

/*
========= Plozia =========
	Author:Plozia
	Problem:P5656 【模板】二元一次不定方程 (exgcd)
	Date:2022/10/10
========= Plozia =========
*/

#include <bits/stdc++.h>
typedef long long LL;

// const int MAXN = ;
LL a, b, c;

int Read()
{
	int sum = 0, fh = 1; char ch = getchar();
	for (; ch < '0' || ch > '9'; ch = getchar()) fh -= (ch == '-') << 1;
	for (; ch >= '0' && ch <= '9'; ch = getchar()) sum = (sum << 3) + (sum << 1) + (ch ^ 48);
	return sum * fh;
}
LL gcd(LL a, LL b) { return (b == 0) ? a : gcd(b, a % b); }
void exgcd(LL a, LL b, LL &x, LL &y)
{
	if (b == 0) { x = 1, y = 0; return ; }
	exgcd(b, a % b, x, y); LL p = x; x = y; y = p - (a / b) * y;
}

void Solve()
{
	a = Read(), b = Read(), c = Read(); LL d = gcd(a, b), x, y;
	if (c % d != 0) { puts("-1"); return ; }
	exgcd(a, b, x, y); x *= c / d; y *= c / d;
	x = (x % (b / d) + (b / d)) % (b / d); if (x == 0) x += b / d; y = (c - a * x) / b;
	if (x > 0 && y <= 0) { printf("%lld %lld\n", x, (((y = (y % (a / d) + (a / d)) % (a / d)) == 0) ? (a / d) : y)); return ; }
	LL Minx = x, Maxy = y; y = (y % (a / d) + (a / d)) % (a / d); if (y == 0) y += a / d; x = (c - b * y) / a;
	LL Maxx = x, Miny = y; printf("%lld %lld %lld %lld %lld\n", (Maxx - Minx) / (b / d) + 1, Minx, Miny, Maxx, Maxy);
}

int main()
{
	int t = Read(); while (t--) Solve(); return 0;
}

3. 总结

exgcd 求 \(ax+by=\gcd(a,b)\) 的方式:类似求 \(\gcd(a,b)\),若 \(bx'+(a\bmod b)y'=\gcd(a,b)\)\((x,y)=(y',x'-\lfloor\dfrac{a}{b}\rfloor y')\)

exgcd 求 \(ax+by=c\) 整数解的一般方式:

  1. \(\gcd(a,b)\not\mid c\) 则无解。
  2. 求出 \(ax+by=\gcd(a,b)\) 的一组特解 \(x',y'\) 后,有 \(x_0=x'\times\dfrac{c}{\gcd(a,b)},y_0=y'\times\dfrac{c}{\gcd(a,b)}\),整数解 \((x,y)\) 的通解如下:

\[\begin{cases}x=x_0+k\times\dfrac{b}{\gcd(a,b)}\\y=y_0-k\times\dfrac{a}{\gcd(a,b)}\end{cases},k\in Z \]

posted @ 2022-04-17 15:05  Plozia  阅读(333)  评论(0编辑  收藏  举报