扩展欧几里德求解线性同余方程组相关~

记录一个菜逼的成长。。

现在理解了，怕以后忘了。

记录一下。。

我们都知道欧几里德算法，又叫辗转相除法。

gcd(a,b) = gcd(b,a%b);

代码实现为

int gcd(int a,int b)
{
return b ? gcd(b,a%b) : a;
}

扩展欧几里德

对于不完全为 0 的非负整数 a，b，gcd（a，b）表示 a，b 的最大公约数，必然

存在整数对 x，y ，使得 gcd（a，b）=ax+by。

代码实现为

int ext_euclid(int a,int b,int &x,int &y)  //求gcd(a,b)=ax+by
{
int t,d;
if (b==0) {x=1;y=0;return a;}
d=ext_euclid(b,a % b,x,y);
//注意交换过程
t=x;//保存下一个方程的x2
x=y;//x是下一个方程的y2，也即是这里的y
y=t-a/b*y;//y = x2 - a / b * y;
return d;
}

求解 x，y的方法的理解

设 a>b。

1，显然当 b=0，gcd（a，b）=a。此时 x=1，y=0；

2，a>b>0 时

设 ax1+ by1= gcd(a,b);

bx2+ (a mod b)y2= gcd(b,a mod b);

根据朴素的欧几里德原理有 gcd(a,b) = gcd(b,a mod b);

则:ax1+ by1= bx2+ (a mod b)y2;

即:ax1+ by1= bx2+ (a - [a / b] * b)y2=ay2+ bx2- [a / b] * by2;

也就是ax1+ by1 == ay2+ b(x2- [a / b] *y2);

根据恒等定理得：x1=y2; y1=x2- [a / b] *y2;//对应函数里的交换过程

这样我们就得到了求解 x1,y1 的方法：x1，y1 的值基于 x2，y2.

上面的思想是以递归定义的，因为 gcd 不断的递归求解一定会有个时候 b=0，所以递归可以结束。

扩展欧几里德求解线性同余方程组

比如 ax = b (modn) 可以表示成 ax + ny = b;

如果方程有解，则满足gcd(a,n)|b ,即 b % gcd(a,n) == 0,

否则无解。

原本应该是这样ax’ + ny’ = gcd(a,n) 因为b是gcd(a,n)的倍数 ,所以可以表示成这样

ax’ * b/gcd(a,n) + ny’ * b/gcd(a,n) = b;(两边同时乘上b/gcd(a,n).)

x = x’ * b / gcd(a,n);

对于同余方程组：

x = a1 (mod m1)

x = a2 (mod m2)

x = a3 (mod m3)

…

x = an (mod mn)

x = m1 * y + a1;

x = m2 * z + a2;

所以可得 m1 * y + a1 = m2 * z + a2;

移项可得: m2 * z - m1 * y = a1 - a2;

即m2*z = a1 - a2 (mod m1);

这里可用扩展欧几里德算法,调用ext_euclid(m2,m1,x’,y’)求解d = gcd(m2,m1);

判断是否有解即(a1 - a2) % gcd(m2,m1) 是否等于0；

如果有解,解为x = x’ * (a1-a2) / d(因为这里是利用m2和a2来求解，所以是用x’来求解);

然后可令a1 = m2 * x + a2(两个同余方程求出来的解)，m1 = m1 * m2 / gcd(m1,m2)(就是最小公倍数)

因为a1可能是负数，所以需要a1 = (a1 % m1 + m1) % m1;

继续跟第三个方程联立求解。

最后的答案就是a1.

中国剩余定理

对于同余方程组：

x = a1 (mod m1)

x = a2 (mod m2)

x = a3 (mod m3)

…

x = an (mod mn)

其中m1、m2…mn互素(如果不互素只能用上面的扩展欧几里德一个一个求并判断是否有解,这里推荐一道题poj2891);

设M = m1 * m2 * … * mn;并设Mi = M / mi,是除了mi以外的n- 1个整数的乘积。

设ti 为Mi mod mi下的逆元 即 ti*Mi = 1 (mod mi)

则通解为 x = a1 * t1 * M1 + a2 * t2 * M2 + … + an * tn * Mn + KM;

在模M意义下，有唯一解，一般我们求的就是这个 x = (a1 * t1 * M1 + a2 * t2 * M2 + … + an * tn * Mn) mod M;

对于求逆元ti 可以调用ext_euclid(mi,Mi,x,y),因为mi 和 Mi是互素的 所以gcd(mi,Mi) = 1.所以求出来的y就是Mi的逆元ti.

代码实现为

typedef long long LL;
LL China(LL W[],LL B[],LL k)   //W为按多少排列，B为剩余个数   W>B  K为组数
{
LL i;
LL d,x,y,a=0,m,n=1;
for (i = 0;i<k;i++)
n *= W[i];
for (i=0;i<k;i++)
{
m=n/W[i];
d=ext_euclid(W[i],m,x,y);
a=(a+y*m*B[i])%n;
}
if (a>0)
return a;
else
return(a+n);
}

PS:还是百度百科看的比较懂。。233