【扩展欧几里得算法】
扩展欧几里得算法（简称扩欧）是欧几里得算法（又称辗转相除法）的扩展。扩欧可以在求得 $a,b$ 的最大公约数的同时，能找到整数 $x,y$，使它们满足裴蜀定理（Bézout's identity，贝祖定理）等式：$ax+by = \text{gcd}(a, b)$。如果 $a$ 是负数，可以把问题转化成 $|a|(-x) + by = \text{gcd}(|a|, b)$，然后令 $x^` = (-x)$。
例如 $a=6,\ b=9$，使用欧几里得算法我们可以得到 $\text{gcd}(6, 9) = 3$，而使用扩欧算法，我们不仅可以计算得到最大公约数 $3$，而且可以得到方程 $6x+9y = 3$ 的整数解 $x,y$。
贝祖定理告诉我们该方程必有整数解，且可以简单证明：
方程 $ax+by = c$ 有整数解  $\leftrightarrow$ $c$ 是 $a$ 和 $b$ 的最大公约数的倍数，即 $c \% \text{gcd}(a, b) \equiv 0$ 成立。
证明：
方程 $ax+by = c$ 有整数解，即 $x,y$ 为整数
$\rightarrow\ ax \% \text{gcd}(a, b) \equiv 0$ 和 $by \% \text{gcd}(a, b) \equiv 0$ 成立
$\rightarrow\ (ax+by) \% \text{gcd}(a, b) \equiv 0$
$\rightarrow\ c \% \text{gcd}(a, b) \equiv 0$
Q.E.D.


对于一个方程 $a∗x+b∗y=\text{gcd}(a,b)$ 来说,我们可以做如下的推导：
设有 $a∗x_1+b∗y_1=\text{gcd}(a,b)$，同时我们有 $b∗x_2+(a\%b)∗y_2=\text{gcd}(b,a\%b)$。
对于这个方程组,我们希望知道的是 $x_1,x_2,y_1,y_2$ 之间的关系，这样我们才可以递归解决这个问题。
我们观察 $b∗x_2+(a\%b)∗y_2$ 这个式子，我们可以将 $(a\%b)$ 写作 $(a−⌊\frac{a}{b}⌋∗b)$，将括号打开常数 $a,b$ 合并，合并之后的结果为 $a∗y_2+b∗(x2−⌊\frac{a}{b}⌋∗y_2))$。
由于欧几里得算法的原理 $\text{gcd}(a,b) \equiv \text{gcd}(b,a\%b)$，我们将两式子联立，对比系数即可得到 $x_1=y_2,y_1=x_2−⌊\frac{a}{b}⌋∗y_2$。
这个递归的边界是什么呢？我们知道，当朴素欧几里得到达边界时，return gcd(a,0)=a，那么边界条件就是对 $a∗x_0+b∗y_0=a$ 求解，很显然此时 $x_0=1,y_0=0$。
当我们递归求出了一个方程的特解时，如何求出这个方程的通解呢？
方程 $a∗x+b∗y=\text{gcd}(a,b)$ 中，如果将 $x$ 加上一个常数 $k_1$，$y$ 减去一个常数 $k_2$，仍然保持原方程成立。那么 $x+k_1,y−k_2$ 就是方程的一个新解。这个 $k$ 应该如何选择呢？
实际上很简单，$a∗(x+k_1)+b∗(y+k_2)=\text{gcd}(a,b)$，打开括号，$a∗x+a∗k_1+b∗y−b∗k_2=\text{gcd}(a,b)$。
保证原方程成立，就需要 $a∗k_1 \equiv b∗k_2$，那么显然 $k_1=b,\ k_2=a$ 是一种合理的情况，但是这样是无法包含所有整数解的，因为我们加上的这个值并非是最小值。
那我们应该加上什么值才行呢？我们发现当 $a∗k_1 \equiv b∗k_2=t∗\text{lcm}(a,b)$ 可以保证得到所有解，于是每次寻找解就可以分别在 $x$ 加上 $\frac{b}{\text{gcd}(a,b)}$，在 $y$ 减去 $\frac{a}{\text{gcd}(a,b)}$，就可以了。
对于方程 $a∗x+b∗y=c$ 我们又该如何求解？
我们发现如果 $(c\%gcd(a,b) \neq 0)$ 那么这个方程是无解的，而如果 $\text{gcd}(a,b)∗t \equiv c$，我们就可以按上面的方法求解之后对我们的解乘上一个 $t\ (t=\frac{c}{\text{gcd}(a,b)})$。




【模板代码】

//求解ax+by=gcd(a,b)
//返回值为 gcd(a,b)
LL exgcd(LL a, LL b, LL &x, LL &y) {
    //出口
    if (b==0) {
        x=1;y=0;
        return a;
    }
    LL gcd=exgcd(b, a%b, y, x);
    y=y-a/b*x;
    return gcd;
}