P1082 [NOIP2012 提高组] 同余方程

问题转化

题目问的是满足 $axmodb=1$ 的最小正整数 $x$。（$a,b$ 是正整数）

但是不能暴力枚举 $x$，会超时。

把问题转化一下。观察 $axmodb=1$，它的实质是 $ax+by=1$：这里 $y$ 是我们新引入的某个整数，并且似乎是个负数才对。这样表示是为了用扩展欧几里得算法。我们将要努力求出一组 $x,y$ 来满足这个等式。稍微再等一下——

问题还需要转化。扩展欧几里得是用来求 $ax+by=gcd(a,b)$ 中的未知数的，怎么牵扯到等于 $1$ 呢？

原理是，方程 $ax+by=m$ 有解的必要条件是 $mmodgcd(a,b)=0$。

---

这个简单证一下。

由最大公因数的定义，可知 $a$ 是 $gcd(a,b)$ 的倍数，且 $b$ 是 $gcd(a,b)$ 的倍数。若 $x,y$ 都是整数，就确定了 $ax+by$ 是 $gcd(a,b)$ 的倍数。

因为 $m=ax+by$，所以 $m$ 必须是 $gcd(a,b)$ 的倍数。

那么 $mmodgcd(a,b)=0$。

---

可得出在这道题中，方程 $ax+by=1$ 的有解的必要条件是 $1modgcd(a,b)=0$。可怜的 $gcd(a,b)$ 只能等于$1$ 了。这实际上就是 $a,b$ 互质。

扩展欧几里得

前提：知道普通欧几里得算法（辗转相除法）。

下面字母挺多，希望你耐心地慢慢地读~

我们拿到了一组 $a,b$。设 $G=gcd(a,b)$。那么目标是求出满足 $ax+by=G$（①）的整数 $x$ 与 $y$。其中 $x$ 应当是满足条件的最小正整数，即答案，而 $y$ 是辅助答案。

（注意，虽然刚刚已经证明本题的 $gcd(a,b)$ 必然等于 $1$，但是扩展欧几里得算法本身过程求的是 $ax+by=gcd(a,b)$ 的解。它正好非常适合本题，不过我们要按照它求解的过程去做）

如果我们先前已经求出了另一组数 $x_2,y_2$，它们满足这么一个式子：$b{x}_2+(amodb)y_2=G$（②），则此时结合①②一定有：$ax+by=b{x}_2+(amodb)y_2$（③）

可见，在这个“如果”实现的时候，我们的目标变成了“求出满足上式的 $x$ 和 $y$”。

其中 $a,b,x_2,y_2$ 都已知，$x,y$ 待求。因为未知数比方程更多，所以没有唯一解。我们先求出一组必然存在的解，最后将在“答案处理”时转为最小解。

怎么求呢？取模运算是 $amodb=a-b\times (a/b)$，所以方程③实际上是：

$ax+by=b{x}_2+(a-b\times (a/b))y_2$

$\Rightarrow ax+by=b{x}_2+a{y}_2-b\times (a/b)y_2$
​
$\Rightarrow ax+by=a{y}_2+b(x_2-(a/b)y_2)$

看上面这个方程，一组必然存在的解出现了：$x=y_2,y=x_2-(a/b)y_2$（④）

可见，我们只要求出 $x_2,y_2$，就能得出正确的 $x,y$。问题是 $x_2,y_2$怎么求。

现在我们手上是 $b,amodb$ 这两个系数，而目标是求出 $x_2$ 和 $y_2$ 满足：$b{x}_2+(amodb)y_2=G$（②）

把①和②对比一下：$ax+by=G$（①）

原方程中的系数 $a$ 变成了②中的系数 $b$，原方程中的系数 $b$ 变成了②中的 $amodb$ 而已。

所以，把新的方程也看作 $ax+by=G$ 的形式（只是系数 $a$ 和 $b$ 的具体数值改变了）。然后按照上面的一模一样下来（其实都只是推导过程），我们发现，最好有 $x_3,y_3$ 来支撑 $x_2,y_2$。

再一模一样下来，我们又需要 $x_4,y_4$ 来支撑 $x_3,y_3$。

……

这个递归中 $a,b$ 不断被替代为 $b,amodb$，这个替换方式与普通欧几里得是一样的，所以最后会出现 $a_n=G,b_n=0$。

这时要直接返回了，我们需要一组 $x_n,y_n$ 满足 $a_n{x}_n+b_n{y}_n=G$（⑤）。

然而该层的 $a_n=G,b_n=0$。所以只要⑤左边取 $x_n=1$，这个方程就妥妥的成立了。

（最后一层的 $y_n$ 建议取 $0$。然而由于 $b=0$，就算返回其它数值，方程也一定成立。但这样的程序容易出错，因为 $y_n$ 在回溯时滚雪球式增长，容易数值越界。下面的代码在最后一层令 $y=7$，开了 long long，通过了此题。）

最后一层结束后，就开始返回，直到最上层。每一层可以轻松地根据下层的 $x_{k+1},y_{k+1}$ 求出当前层的 $x_k,y_k$。

整个过程就是：以辗转相除的方式向下递进，不断缩小系数，保证会出现有确定解的最后一层。

答案处理

一个遗留问题：我们将要求出来的 $x,y$ 仅仅保证满足 $ax+by=1$，而 $x$ 不一定是最小正整数解。有可能太大，也有可能是负数。这依然可以解决，那就是，$x$ 批量地减去或加上 $b$，能保证 $ax+by=1$，因为：

$ax+by=1$

$ax+by+k\times ba-k\times ba=1$

$a(x+kb)+(y-ka)b=1$

1.显然这并不会把方程中任何整数变成非整数。

2.“加上或减去 $b$”不会使 $x$ 错过任何解。可以这么理解：

---

已经求出一组整数 $x,y$ 使得 $ax+by=1$，也就是 $\frac{1-ax}{b}=y$。$y$ 是整数，可见目前 $1-ax$ 是 $b$ 的倍数。

现在想改变 $x$ 并使得方程仍然成立。已知 $a,b$ 互质，假若 $x$ 的变化量（$\Delta x$）不是 $b$ 的倍数，则 $1-ax$ 的变化量（$-a\times \Delta x$）也不是 $b$ 的倍数，这会使得 $1-ax$ 不再是 $b$ 的倍数，则 $y$ 不是整数了。

仅当 $x$ 的变化量是 $b$ 的倍数时，$1-ax$ 能保持自己是 $b$ 的倍数，此时就出现新的解了。

因此到最后，如果 $x$ 太小就不断加 $b$ 直到大于等于 $0$，太大则一直减 $b$，直到最小正整数解。也就是这么写：

x = (x % b + b) % b;//括号中取模再加，可以处理负数

代码

推导中的所有 $x,y$ 共用全局变量 long long x, y，传递也很方便。

#include<bits/stdc++.h>
using namespace std;
long long x, y;//目前方程真正的解 
void exgcd(long long a, long long b)
{
    //当前目的：求解 ax + by = gcd(a, b) 这么一个方程
    if(b == 0) //a, b不断改变的过程中，b最终必然会成为0
    {
        //在 b = 0 时方程还要成立？ 使 x = 1, y = 0 ，必然成立 
        x = 1;
        y = 7; //建议返回0。不过y = 7能AC，证明了最后一个等式不受最后一个y影响
        return;
    } 
    exgcd(b, a % b);//把下一层系数传进去（先求下一个方程的解 ）
    //现在我们已经拿到了下一个方程的解x, y
    long long tx = x;//暂时存一下x，别丢了
    x = y;
    y = tx - a / b * y; 
}
int main()
{
    long long a, b;
    cin >> a >> b;
    exgcd(a, b);
    x = (x % b + b) % b;//我们求出来的x必然满足方程，但不一定是最小正整数解，所以要进行答案处理
    printf("%lld\n", x);
    return 0;
}

求这样一个满足 $axmodb=1$ 的 $x$ 有什么用呢？一个重要用途如下。

这个 $x$ 就是：$a$ 在模 $b$ 意义下的乘法逆元。

在模 $b$ 意义下，如果想要除以 $a$ 就非常麻烦。这时候乘以 $a$ 的逆元等效于除以 $a$。

假设我们已经算出了：

ans = (a * b * c) % p;

现在要从 $ans$ 中把 $b$ 给除掉，如何处理 $ans$？

如果像下面这样直接除会出错（举个实例会很直观）：

ans = ans / b % p;

所以我们就求出 $b$ 的逆元 $x$（满足 $bxmodp=1$），然后直接算这个：

ans = ans * x % p;

原理可以这么理解：$abcx\equiv ac(bx)\equiv ac(modp)$