【算法讲2：拓展欧几里得（简略讲）】求解 ax+by=c

拓展欧几里得算法

欧几里得算法简介
拓展欧几里得算法简介
- 【递归方法：倒序带入法】
- 【递归代码】
- 【递推方法：滚动变量记录法】
- 【递推代码】
应用：如何解 ax+by=c
- 测试代码

欧几里得算法简介

大部分参考与查阅：《初等数论及其基本应用》原书第六版。
指辗转相除法。
令整数 r0=a，r1=br_0=a，r_1=br0=a，r1=b，满足 a≥b>0a\ge b>0a≥b>0，连续做除法得到：
rj=rj+1qj+1+rj+2（j=0ton−2），rn+1=0r_j=r_{j+1}q_{j+1}+r_{j+2}（j=0\ to\ n-2），r_{n+1}=0rj=rj+1qj+1+rj+2（j=0 to n−2），rn+1=0
那么（a,b）= rnr_nrn
例如对于 (252,198)的例子：

j	rjr_jrj	rj+1r_{j+1}rj+1	qj+1q_{j+1}qj+1	rj+2r_{j+2}rj+2
0	252	198	1	54
1	198	54	3	36
2	54	36	1	18
3	36	18	2	0

算法的时间复杂度（拉梅定律）： O(5log⁡10b)O(5\log_{10} b)O(5log10b)

拓展欧几里得算法简介

在欧几里得算法的基础上，我们记录出一些数值为了
用线性组合来表示最大公因子
比如，(252,198)=18(252,198)=18(252,198)=18，但是我想知道 252s+198t=18252s+198t=18252s+198t=18，其中两个系数该怎么取。

【递归方法：倒序带入法】

见上倒数第二步，18=54−1∗3618=54-1*3618=54−1∗36
前一步为 36=198−3∗5436=198-3*5436=198−3∗54
带入得到 18=54−1∗(198−3∗54)=4∗54−1∗19818=54-1*(198-3*54)=4*54-1*19818=54−1∗(198−3∗54)=4∗54−1∗198
⋯\cdots⋯

【算法描述】：
原理：更新 x、yx、yx、y，使得 rj∗x+rj+1∗y=(a,b)r_j*x+r_{j+1}*y=(a,b)rj∗x+rj+1∗y=(a,b)恒成立。

步骤：首先跑一边欧几里得算法，接下来逆序计算。
首先，rn=(a,b),rn+1=0r_n=(a,b),r_{n+1}=0rn=(a,b),rn+1=0
对于特定的 jjj ，如果已经求得 (a,b)=srj+trj−1(a,b)=sr_j+tr_{j-1}(a,b)=srj+trj−1
那么，因为 rj=rj+1qj+1+rj+2（j=0ton−2），rn+1=0r_j=r_{j+1}q_{j+1}+r_{j+2}（j=0\ to\ n-2），r_{n+1}=0rj=rj+1qj+1+rj+2（j=0 to n−2），rn+1=0
等式变一变，我们有：
rj=rj−2−rj−1qj−1r_j=r_{j-2}-r_{j-1}q_{j-1}rj=rj−2−rj−1qj−1
带入得到：
(a,b)=s(rj−2−rj−1qj−1)+trj−1=(t−sqj−1)rj−1+srj−2\begin{aligned}(a,b)&=s(r_{j-2}-r_{j-1}q_{j-1})+tr_{j-1}\\&=(t-sq_{j-1})r_{j-1}+sr_{j-2}\end{aligned}(a,b)=s(rj−2−rj−1qj−1)+trj−1=(t−sqj−1)rj−1+srj−2
如果仔细看一下式子，便可以发现：

经过一次变换之后，
sss变成(t−sqj−1)(t-sq_{j-1})(t−sqj−1)
ttt变成sss

经过多次变换之后，我们的 s、rs、rs、r的系数就可以从rn、rn+1r_n、r_{n+1}rn、rn+1 变成 r0、r1r_0、r_1r0、r1，也就是我们想要得到的线性组合的两个系数了。

【递归代码】

ll ex_gcd(ll a,ll b,ll& x,ll& y) { if(b==0) {x=1;y=0;return a;}ll g=ex_gcd(b,a % b,x,y);ll tmp = x;x = y; y = tmp - a / b * y;return g;
}

【递推方法：滚动变量记录法】

原理：我们保持 rj=sja+tjbr_j=s_ja+t_jbrj=sja+tjb 其恒成立

首先，我们递归定义一下：

(a,b)=sna+tnb(a,b)=s_na+t_nb(a,b)=sna+tnb

s0=1，t0=0s_0=1，t_0=0s0=1，t0=0
s1=0，t1=1s_1=0，t_1=1s1=0，t1=1

∀j≥2sj=sj−2−qj−1sj−1tj=tj−2−qj−1tj−1\forall\ j\ge 2\\s_j=s_{j-2}-q_{j-1}s_{j-1}\\t_j=t_{j-2}-q_{j-1}t_{j-1}∀ j≥2sj=sj−2−qj−1sj−1tj=tj−2−qj−1tj−1

证明（第二数归，略）

【递推代码】

ll ex_gcd2(ll a,ll b,ll& s1,ll &t1) {s1 = 0;t1 = 1;ll s2 = 1,t2 = 0;ll q;bool fir = true;while(b){if(!fir){ll n1 = s2 - q * s1;ll n2 = t2 - q * t1;s2 = s1;s1 = n1;t2 = t1;t1 = n2;}if(b)q = a / b;ll tmp = a;a = b;b = tmp % b;fir = false;}return a;
}

应用：如何解 ax+by=c

【题目介绍】
a、b、ca、b、ca、b、c是已给出的整数。
我想求得一组整数解 x,yx,yx,y，满足 a∗x+b∗y=ca*x+b*y=ca∗x+b∗y=c

首先，我们得到 g=(a,b)g=(a,b)g=(a,b)，并且容易知道 g∣(ax+by)g|(ax+by)g∣(ax+by)
换句话说，只有 g∣cg|cg∣c ，该方程才有解。
并且，若 g∣cg|cg∣c，该方程有无穷多个解。

【常规解法】：

求出 g=(a,b)g=(a,b)g=(a,b)
若 g∤cg\nmid cg∤c，则该方程无解，直接跳出。否则有无穷多个解。
我们先求出 ax+by=gax+by=gax+by=g的解，然后方程两边同时乘以 cg，即可得到解\frac{c}{g}，即可得到解gc，即可得到解
解方程 ax+by=gax+by=gax+by=g ，首先得求得该方程的一个解，即求出 (a,b)(a,b)(a,b) 的线性组合，拓展欧几里得求出的两个系数即为其方程的一个解 x0,y0x_0,y_0x0,y0
方程 ax+by=gax+by=gax+by=g 的通解为：
{x=x0+kbgy=y0−kagk∈Z\begin{cases} x=x_0+k\frac{b}{g}\\ y=y_0-k\frac{a}{g}\\ \end{cases} \\k\in\mathbb{Z} {x=x0+kgby=y0−kgak∈Z
方程 ax+by=cax+by=cax+by=c 的通解为：
{x=x0cg+kbgy=y0cg−kagk∈Z\begin{cases} x=x_0\frac{c}{g}+k\frac{b}{g}\\ y=y_0\frac{c}{g}-k\frac{a}{g}\\ \end{cases} \\k\in\mathbb{Z} {x=x0gc+kgby=y0gc−kgak∈Z

测试代码

可以看递归方法是否满足 rj∗x+rj+1∗y=(a,b)r_j*x+r_{j+1}*y=(a,b)rj∗x+rj+1∗y=(a,b)恒成立。
可以看递推方法是否满足 rj=sja+tjbr_j=s_ja+t_jbrj=sja+tjb 其恒成立
可以看是否对于任意的整数 kkk 都满足 ax+by=cax+by=cax+by=c
只要在宏定义上面修改A、B、C的值即可。（可以看看 C∤(A,B)C\nmid(A,B)C∤(A,B)会发生什么）

/*_            __   __          _          _
| |           \ \ / /         | |        (_)
| |__  _   _   \ V /__ _ _ __ | |     ___ _
| '_ \| | | |   \ // _` | '_ \| |    / _ \ |
| |_) | |_| |   | | (_| | | | | |___|  __/ |
|_.__/ \__, |   \_/\__,_|_| |_\_____/\___|_|__/ ||___/
*/
#define ll long long
ll A = 252;
ll B = 198;
ll C = 36;
#define show(x) std::cerr << #x << " = " << x << std::endl
#define show2(x,y) std::cerr << #x << " = " << x << "   " << #y << " = " << y << std::endl
#define show3(x,y,z) std::cerr << #x << " = " << x << "   " << #y << " = " << y << "   " << #z << " = " << z << std::endl
#define show4(x,y,z,p) std::cerr << #x << " = " << x << "   " << #y << " = " << y << "   " << #z << " = " << z << "   " << #p << " = " << p << std::endl
ll ex_gcd(ll a,ll b,ll& x,ll& y) { if(b==0) { x=1;y=0;show2(x,y); return a; } ll g=ex_gcd(b,a%b,x,y); ll tmp=x; x=y; y=tmp-a/b*y; show3(x,y,a * x + b * y);return g; }ll ex_gcd2(ll a,ll b,ll& s1,ll &t1) {s1 = 0;t1 = 1;ll s2 = 1,t2 = 0;ll q;bool fir = true;while(b){if(!fir){ll n1 = s2 - q * s1;ll n2 = t2 - q * t1;s2 = s1;s1 = n1;t2 = t1;t1 = n2;show2(s1,t1);show(s1 * A + t1 * B);}if(b)q = a / b;ll tmp = a;a = b;b = tmp % b;fir = false;}return a;
}int main()
{ll x,y;cout << ex_gcd(A,B,x,y) << endl;show2(x,y);puts("----------");ll a,b;cout << ex_gcd2(A,B,a,b) << endl;show2(x,y);ll g = __gcd(A,B);for(int k = -5;k <= 5;++k){ll xx = x * C / g + k * B / g;ll yy = y * C / g - k * A / g;show2(xx,yy);cout << xx * A + yy * B << endl;}return 0;
}