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

欧几里得算法简介

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

j	$r_j$	$r_{j+1}$	$q_{j+1}$	$r_{j+2}$
0	252	198	1	54
1	198	54	3	36
2	54	36	1	18
3	36	18	2	0

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

拓展欧几里得算法简介

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

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

见上倒数第二步，$18=54-1\ast 36$
前一步为 $36=198-3\ast 54$
带入得到 $18=54-1\ast (198-3\ast 54)=4\ast 54-1\ast 198$
$\cdots$

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

步骤：首先跑一边欧几里得算法，接下来逆序计算。
首先，$r_n=(a,b),r_{n+1}=0$
对于特定的 $j$ ，如果已经求得 $(a,b)=s{r}_j+t{r}_{j-1}$
那么，因为 $r_j=r_{j+1}{q}_{j+1}+r_{j+2}（j=0ton-2），r_{n+1}=0$
等式变一变，我们有：
$r_j=r_{j-2}-r_{j-1}{q}_{j-1}$
带入得到：
$\begin{array}{rl}(a,b)& =s(r_{j-2}-r_{j-1}{q}_{j-1})+t{r}_{j-1}\\ & =(t-s{q}_{j-1})r_{j-1}+s{r}_{j-2}\end{array}$
如果仔细看一下式子，便可以发现：

经过一次变换之后，
$s$变成$(t-s{q}_{j-1})$
$t$变成$s$

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

【递归代码】

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; 
}

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

原理：我们保持 $r_j=s_{j}a+t_{j}b$ 其恒成立

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

$(a,b)=s_{n}a+t_{n}b$

$s_0=1，t_0=0$
$s_1=0，t_1=1$

$$\begin{gathered} \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} \end{gathered}$$

证明（第二数归，略）

【递推代码】

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、c$是已给出的整数。
我想求得一组整数解 $x,y$，满足 $a\ast x+b\ast y=c$

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

【常规解法】：

1. 求出 $g=(a,b)$
2. 若 $g\nmid c$，则该方程无解，直接跳出。否则有无穷多个解。
3. 我们先求出 $ax+by=g$的解，然后方程两边同时乘以 $\frac{c}{g}，即可得到解$
4. 解方程 $ax+by=g$ ，首先得求得该方程的一个解，即求出 $(a,b)$ 的线性组合，拓展欧几里得求出的两个系数即为其方程的一个解 $x_0,y_0$
5. 方程 $ax+by=g$ 的通解为 ：
   $$\begin{gathered} \{\begin{array}{l}x=x_0+k\frac{b}{g}\\ y=y_0-k\frac{a}{g}\end{array} \\ k\in \mathbb{Z} \end{gathered}$$
6. 方程 $ax+by=c$ 的通解为 ：
   $$\begin{gathered} \{\begin{array}{l}x=x_0\frac{c}{g}+k\frac{b}{g}\\ y=y_0\frac{c}{g}-k\frac{a}{g}\end{array} \\ k\in \mathbb{Z} \end{gathered}$$

测试代码

可以看递归方法是否满足 $r_j\ast x+r_{j+1}\ast y=(a,b)$恒成立。
可以看递推方法是否满足 $r_j=s_{j}a+t_{j}b$ 其恒成立
可以看是否对于任意的整数 $k$ 都满足 $ax+by=c$
只要在宏定义上面修改A、B、C的值即可。（可以看看 $C\nmid (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;
}