P1082 [NOIP2012 提高组] 同余方程
问题转化
题目问的是满足 \(ax \bmod b = 1\) 的最小正整数 \(x\)。(\(a,b\) 是正整数)
但是不能暴力枚举 \(x\),会超时。
把问题转化一下。观察 \(ax \bmod b = 1\),它的实质是 \(ax+by=1\):这里 \(y\) 是我们新引入的某个整数,并且似乎是个负数才对。这样表示是为了用扩展欧几里得算法。我们将要努力求出一组 \(x,y\) 来满足这个等式。稍微再等一下——
问题还需要转化。扩展欧几里得是用来求 \(ax+by=gcd(a,b)\) 中的未知数的,怎么牵扯到等于 \(1\) 呢?
原理是,方程 \(ax+by=m\) 有解的必要条件是 \(m \bmod gcd(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)\) 的倍数。
那么 \(m \bmod gcd(a,b)=0\)。
可得出在这道题中,方程 \(ax+by=1\) 的有解的必要条件是 \(1 \bmod gcd(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\),它们满足这么一个式子:\(bx_2+(a \bmod b)y_2=G\)(②),则此时结合①②一定有:\(ax+by=bx_2+(a \bmod b)y_2\)(③)
可见,在这个“如果”实现的时候,我们的目标变成了“求出满足上式的 \(x\) 和 \(y\)”。
其中 \(a,b,x_2,y_2\) 都已知,\(x,y\) 待求。因为未知数比方程更多,所以没有唯一解。我们先求出一组必然存在的解,最后将在“答案处理”时转为最小解。
怎么求呢?取模运算是 \(a \bmod b=a-b \times (a/b)\),所以方程③实际上是:
\(ax+by=bx_2+(a-b \times (a/b))y_2\)
\(\Rightarrow ax+by=bx_2+ay_2-b \times (a/b)y_2\)
\(\Rightarrow ax+by=ay_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, a \bmod b\) 这两个系数,而目标是求出 \(x_2\) 和 \(y_2\) 满足:\(bx_2+(a \bmod b)y_2=G\)(②)
把①和②对比一下:\(ax+by=G\)(①)
原方程中的系数 \(a\) 变成了②中的系数 \(b\),原方程中的系数 \(b\) 变成了②中的 \(a \bmod b\) 而已。
所以,把新的方程也看作 \(ax+by=G\) 的形式(只是系数 \(a\) 和 \(b\) 的具体数值改变了)。然后按照上面的一模一样下来(其实都只是推导过程),我们发现,最好有 \(x_3,y_3\) 来支撑 \(x_2,y_2\)。
再一模一样下来,我们又需要 \(x_4,y_4\) 来支撑 \(x_3,y_3\)。
……
这个递归中 \(a,b\) 不断被替代为 \(b, a \bmod b\),这个替换方式与普通欧几里得是一样的,所以最后会出现 \(a_n=G, b_n=0\)。
这时要直接返回了,我们需要一组 \(x_n, y_n\) 满足 \(a_nx_n+b_ny_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×ba−k×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\) 的变化量(\(Δx\))不是 \(b\) 的倍数,则 \(1-ax\) 的变化量(\(−a \times Δ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;
}
求这样一个满足 \(ax \bmod b=1\) 的 \(x\) 有什么用呢?一个重要用途如下。
这个 \(x\) 就是:\(a\) 在模 \(b\) 意义下的乘法逆元。
在模 \(b\) 意义下,如果想要除以 \(a\) 就非常麻烦。这时候乘以 \(a\) 的逆元等效于除以 \(a\)。
假设我们已经算出了:
ans = (a * b * c) % p;
现在要从 \(ans\) 中把 \(b\) 给除掉,如何处理 \(ans\)?
如果像下面这样直接除会出错(举个实例会很直观):
ans = ans / b % p;
所以我们就求出 \(b\) 的逆元 \(x\)(满足 \(bx \bmod p=1\)),然后直接算这个:
ans = ans * x % p;
原理可以这么理解:\(abcx≡ac(bx)≡ac(\bmod p)\)