TAT,感觉这是2012年NOIP唯一的一道我可以做的题。12年题都太神不忍直视,先写道简单的安慰一下自己。
ax ≡ 1 (modb),其实x就是a的逆元。可以用扩展欧几里得算法来求。我知道的还有一种方法用欧拉函数,应该也可以用来接这道题。不过那个方法还需要筛素数,求欧拉函数,如果像这道题一样仅仅求个逆元是比扩展欧几里得算法慢的。
然而在此之前我还不会写扩展欧几里得,算是借这道题学习了一下。
很多人都知道欧几里得算法,用来求最大公约数gcd(a,b)。扩展欧几里得是用来解一类方程:ax+by+c=0的整数解。转化为ax+by = gcd(a,b),利用数学归纳法解此方程。可知,当-c是gcd(a,b)的整倍数的时候原方程有整数解。
以下转载自(详见算法导论):
http://www.cppblog.com/mythit/archive/2009/06/12/87514.html
推论1:方程ax=b(mod n)对于未知量x有解,当且仅当gcd(a,n) | b。
推论2:方程ax=b(mod n)或者对模n有d个不同的解,其中d=gcd(a,n),或者无解。
定理1:设d=gcd(a,n),假定对整数x和y满足d=ax+by(比如用扩展Euclid算法求出的一组解)。如果d | b,则方程ax=b(mod n)有一个解x0满足x0=x*(b/d) mod n 。特别的设e=x0+n,方程ax=b(mod n)的最小整数解x1=e mod (n/d),最大整数解x2=x1+(d-1)*(n/d)。
定理2:假设方程ax=b(mod n)有解,且x0是方程的任意一个解,则该方程对模n恰有d个不同的解(d=gcd(a,n)),分别为:xi=x0+i*(n/d) mod n 。
接下来就要膜拜一下刘汝佳叔叔,C++的Exgcd实现只用2行:
void Exgcd(int a, int b, int &d, int &x, int &y){
if(!b) d = a, x = 1, y = 0;
else { Exgcd(b, a%b, d, y, x); y -= a/b*x;}
}
不过方程的解是不唯一的,我们如何由函数递归解出的x0,y0得到xn,yn呢?
a*x0+b*y0 = a*x1+b*y1;
a*(x0-x1) = b*(y1-y0);
a = a/gcd(a,b), b = b/gcd(a,b);
a*(x0-x1) = b*(y1-y0);
因为此时a,b互质,所以(x0-x1)必然是b的整倍数,(y1-y0)是a的整倍数。即:x0 = k*b+x1, y1 = k*a+y0(k为整数)。由此式可以得到其他解。
下面要做的就是从ax ≡ 1 (mod b)推导到ax+by+c=0这种形式:
ax % b = 1
ax = k*b + 1
ax + b*y = 1 (y == -k)
题目又说必然有解,所以gcd(a, b)必然等于1。
#include <cstdio>
void Exgcd(int a, int b, int &d, int &x, int &y){
if(!b) d = a, x = 1, y = 0;
else {Exgcd(b, a%b, d, y, x); y -= a/b*x;}
}
int main()
{
int a, b, d, x, y;
scanf("%d %d", &a, &b);
Exgcd(a, b, d, x, y);
while(x > 0) x -= b;
while(x < 0) x += b; //gcd(a,b)==1 ∴b == b/gcd(a,b)
printf("%d", x);
return 0;
}