【NOIP2012提高组】同余方程 （扩展欧几里得算法）

【NOIP2012提高组】同余方程 (Standard IO)

Time Limits: 1000
ms  Memory Limits: 131072
KB  Detailed Limits  

Description

求关于x的同余方程ax ≡ 1 (mod b)的最小正整数解。 

Input

输入文件为mod.in。 

输入只有一行，包含两个正整数 a, b，用一个空格隔开。 

Output

输出文件为mod.out。 

输出只有一行，包含一个正整数 x ，即最小正整数解。输入数据保
证一定有解。 

Sample Input

3 10

Sample Output

7

Data Constraint

Hint

对于40%的数据，2 ≤b≤ 1,000； 

对于60%的数据，2 ≤b≤ 50,000,000； 

对于100%的数据，2 ≤a, b≤ 2,000,000,000。 

扩展欧几里得算法。

扩展欧几里得算法是指对于两个数a,b。一定能找到x,y（均为整数，但不满足一定是正数）满足x*a+y*b=gcd(a,b).gcd（x,y)是指x 与 y的最大公约数。

也就是说，找到了一个特解，所有解都可以求出。

怎么求呢？

假设x1 与 y1就是一组特解，那么全部解就有：

我们知道，当被减数与减数同时加上一个数，那么他们的差是不变的。(x1*a+y1*b=x1*a-(-y1*b)=gcd(x,y))

现在，我们要保证xn*a与—yn*b的差也要是gcd(a,b)，那么我们得对a*x1与-y1*b同时加上一个数，又要保证加完这个数后，他还是a和b的倍数，就是说，这个数要是a的倍数，也要是b的倍数。那么加上的数只能是a与b的最小公倍数的倍数。

总结一下，

就是说：

所有的解为：

Xn=X1+b*t/gcd(a,b).

Yn=-Y1+a*t/gcd(a,b).

为什么是这些算式，请自行思考。

现在的问题转移到求他的特解了。我们发现当求辗转相除法到底的时候，b=0，那么这时的特解就为1与0了。这时，考虑上一层递归的特解情况。

看看辗转相除法：

R(a,b)=R(b,a mod b);

假设当前特解为X，Y。而他的下一层递归的特解为X1,Y1。X1与Y1我们是已经知道的了。

那么就有：

X*a+Y*b=X1*b+Y1*(a mod b)

a mod b=a-a div b*b.

右边的式子=X1*b+Y1*(a-a div b*b)

                 =X1*b+Y1*a-Y1*（a div b)*b

 =Y1*a+b*(X1-Y1*(a div b))

X*a+Y*b=Y1*a+（X1-Y1*(a div b))*b。

转的过程为：

X=Y1;

Y=X1-
a273
Y1*(a div b)。

那么就可以求导出一个特解了（当递归结束时）。

然后就可以进行求解了。

上题也如此。数据一定保证a与b的最大公约数一定是1.不然无解。

他是要求(a*x) mod b=1.那么可以推出公式：a*x=b*y+1.也就是： a*x-b*y=1。好像有点奇怪。-？其实我们可以设C=-y,那么问题圆满解决。变：a*x+c*y=1.

所以求出特解A与C。

答案必须是正整数。其实也没有什么难度，可以套用上边的求Xn的通用公式，也可以用这种方法：

求出特解后用MOD法求出最小整数。算出特解的a*x，判断他是不是正数，如果是MOD（a*b)，否则则+a*b*((-a*x) div (a*b)+1).原因请自行思考。最后别忘了DIV A

代码如下：

var

        a,b,x,y,sum,ans:int64;

procedure find(a,b:longint);

var

        tx,ty:longint;

begin

        if b=0 then

          begin

            x:=1;

            y:=0;

            exit;

          end

        else

          find(b,a mod b);

        tx:=x;

        ty:=y;

        x:=ty;

        y:=tx-a div b*ty;

end;

begin

        assign(input,'mod.in');

        assign(output,'mod.out');

        reset(input);

        rewrite(output);

        read(a,b);

        find(a,b);

        sum:=x*a;

        if sum>0 then

          ans:=sum mod (a*b) div a

        else

          begin

            sum:=sum+a*b*((0-sum) div (a*b)+1);

            ans:=sum div a;

          end;

        writeln(ans);

        close(input);

        close(output);

end.