欧几里得/扩展欧几里得/模线性方程组(中国剩余定理 以及mod不互质的情况)
/***********欧几里得算法**************/
辗转相除法求最大公约数。GCD(a,b)=GCD(b,a mod b)稍微证明一下:
(参考:算法导论)证明的思路是大致是这样的:
证明 GCD(a,b) | GCD(b,a mod b) 并且 GCD(b,a mod b) | GCD(a,b)
先证GCD(a,b) | GCD(b,a mod b):设 d=GCD(a,b),则d|a 并且 d|b那么
a=r+kb (k为整数)也就是说 (a mod b) 是a和b的线性组合,
所以有d|(a mod b)因为 d|b 并且 d|(a mod b)
所以 d|GCD(b,a mod b)也就是 GCD(a,b) | GCD(b,a mod b)类似的,GCD(b,a mod b) | GCD(a,b)
//递归 int gcd(int a,int b)
{
return b?gcd(b,a%b):a;
}
//非递归 int gcd(int a,int b)
{
int t;
while(b)
{
t=a;a=b;b=t%b;
}
return a;
}
/**********扩展欧几里得算法*************/
用这个算法可以求 X mod b=m 这种模线性方程(求的是X的通解,整数解)。
变形一下,得 X+by=m,其中y为整数。通常写成 ax+by=m 的形式,然后解得x的通解,进而得到X=ax的通解。
其实就是在求GCD(a,b)的同时把解出满足 ax+by=m 的(x,y)通解。首先,a和b的线性组合能表示的最
小的正整数为GCD(a,b),也就是说ax+by=m,m的最小值为GCD(a,b)ax+by=m有整数解的前提是 GCD(a,b)|m
设d=GCD(a,b)用EXGCD(待会儿再来讲这个函数)
可以解出ax+by=d的一个整数解(x‘,y’),则(x‘*m/d,y’*m/d)就是ax+by=m的一个整数解。
令(x0,y0)=(x‘*m/d,y’*m/d)那么,有(x0+b/d*t,y0-a/d*t)是ax+by=m的通解。
证一下:
由已知得ax0+by0=m将(x,y)=(x0+b/d*t,y0-a/d*t)代入ax+by,
化简一下就有ax0+ab/d*t+by0-ab/d*t=ax0+by0=m证毕。
所以只要得到ax+by=m的一组解,就可以用(x,y)=(x0+b/d*t,y0-a/d*t)得到满足要求的解。
一般要求的是满足条件的最小正整数x。
可以用x=(x0 mod (b/d) + (b/d) ) mod (b/d)
如果求的是X mod b=m的最小正整数解,则X=(ax mod b + b) mod b。这么写是为了避免a为负数的情况下造成的错误。
然后,大家肯定很纳闷这个(x',y')到底怎么求,接下来就讲的是用EXGCD 解 ax+by=d的一个特解的求法。
原来说过的,GCD(a,b)=GCD(b,a mod b)。
那么有ax+by=d………………………………………………(1)
a'x'+b'y'=d……………………………………(2)
其中(a',b')=(b,a mod b)可设 a=bk+t……………………………………(3)
将(3)代入(1),有(bk+t)x+by=d整理,
得b(kx+y)+tx=d………………………………(4)
对比(2)和(4),且用(a',b')替换(b,t),
有a'x' + b'y' =da'(kx+y) + b'x =d对比两式的系数,
有(x',y')=(kx+y,x)整理下,有(x,y)=(y',x'-ky')
也就是说,当深层的(x',y')解出后,就能推出外层的(x,y)边界条件是 b=0时 d=a
那么 ax+by=d 的解(x,y)=(1,0)用EXGCD就能递归地解得 ax+by=d 的一个特解 (x',y')
int Extended_Euclid(int a,int b,int &x,int &y)
{
if(!b)
{x=1;y=0;return a;}
int d = Extended_Euclid(b,a%b,x,y);
int t = x;
x=y;
y=t-(a/b)*y;
return d;
}
二元一次不定方程求解:ax+by=c
gcd(a,b)|c是有整数解的充分必要条件
通过欧几里德扩展定理Extende_Euclid(a,b,x,y)求出一组特殊解:x0=x*c/gcd(a,b)、y0=y*c/gcd(a,b)
通解X=x0+t*b/gcd(a,b)、Y=y0-t*a/gcd(a,b);(t为正整数)
MODULAR-LINEAR_EQUATION_SOLVER(a, b, n)
(d,x’,y’) ← EXTENDED-EUCLID(a, n)
if (d | b)
then x0 ← x’(b/d)mod n
for i ← 0 to d-1
do print(x0 + i(n / d)) mod n
else print “no solution”
求解模线性方程:ax=b (mod n)
通过欧几里德扩展定理Extended_Euclid(a,n,x,y)求出一个特解:x0=x*b/gcd(a,n)%n;
通解X=x0+i*n/gcd(a,n);
/**************线性同余方程**************/
问题是,求满足以下模线性方程组的 X(通常求的是X的最小正整数解)
X mod m1=r1
X mod m2=r2
...
...
...
X mod mn=rn其中,mi之间两两互质线性同余方程中特殊的(用中国剩余定理)
m=m1*m2*.......*mn;
Mi=m/mi;
X=M1’M1*b1+M2’M2*b2+.......+Mn’*Mn*bn;
其中Mi’为乘率:
Mi’*Mi=1(mod mi);
可以通过Extended_Euclid(mi,M,x,y)得出Mi’=y;
int Extended_Euclid(int a, int b, int & x, int & y)
{
if (b==0)
{ x=1; y=0; return a; }
int d = Extended_Euclid(b, a%b, x, y);
int t = x; x = y; y = t - a / b * y;
return d;
}
int china( int n)
{
int i, d, x, y, M, sum, n = 1;
Sum=0;
for (i = 0; i < n; i++) n *= m[i];
for (i = 0; i < n; i++)
{
M= n / m[i];
d = Extended_Euclid(m[i], M, x, y);
sum = (sum+ y * m * b[i]) % n;
}
if (sum > 0) return sum;
else return (sum + n);
}
用中国剩余定理有一个前提条件,就是mi之间是两两互质的!如果mi并不满足两两互质的话。。怎么解呢?
/*****************一般模线性方程*********************/
同样是求这个东西。。
X mod m1=r1
X mod m2=r2
...
...
...
X mod mn=rn
首先,我们看两个式子的情况
X mod m1=r1……………………………………………………………(1)
X mod m2=r2……………………………………………………………(2)
则有 X=m1*k1+r1………………………………………………………………(*)
X=m2*k2+r2那么 m1*k1+r1=m2*k2+r2整理,
得m1*k1-m2*k2=r2-r1令(a,b,x,y,m)=(m1,m2,k1,k2,r2-r1),原式变成ax+by=m熟悉吧?
解线性同余方程组(中国剩余定理中的两两不互质的情况)
int Extended_Euclid(int a,int b,int &x,int &y)
{
if(!b)
{x=1;y=0;return a;}
int d = Extended_Euclid(b,a%b,x,y);
int t = x;
x=y;
y=t-(a/b)*y;
return d;
}
int CRT(int n,int &LCM)
{
int i,x,y;
int d,c,t,mm=m[0],bb=b[0];
for(i=1;i<M;i++)
{
d=exgcd(mm,m[i],x,y);
c=b[i]-bb;
if(c%d)
return -1;
t=m[i]/d;
bb=mm*((c/d*x%t+t)%t)+bb;
mm=mm*m[i]/d;
}
LCM=mm;
return bb;
}