模线性方程组（转载）

/**********************欧几里得算法***********************/

其实很早就知道这个东西了。。也就是以前所说的辗转相除法求最大公约数。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)

大家都会证这个了吧？真聪明。。。。。。。。。

代码我试了三个，其中第三个是算法导论那章末尾的习题。看着觉得很NB，不过稍微测试了下（10^5个int，10^6个int，10^6个long long）。。发现效率也优不了多少，时间都差不多。。也许是我没写好吧。。
建议一般的可以直接用递归的版本（第一种）。非递归版本（第二种）是以防万一。至于位运算的那个版本（第三种）。。。。以后再查查看怎么写才有效率吧。。

//Least Common Multiple
顺便一提，LCM(a,b)=a*b/GCD(a,b)。这在求解模线性方程组时有一定的用处。
这个可以用唯一分解定理证。

a和b都可唯一分解成以下形式：
a=p1^s1*p2^s2*...*pn^sn
b=p1^t1*p2^t2*...*pn^tn
其中，pi为质数，si和ti为非负整数。那么
GCD(a,b)=p1^min(s1,t1)*p2^min(s2,t2)*...*pn^min(sn,tn)
LCM(a,b)=p1^MAX(s1,t1)*p2^max(s2,t2)*...*pn^max(sn,tn)
显然有
GCD(a,b)*LCM(a,b)=a*b
所以，LCM(a,b)=a*b/GCD(a,b)
证毕..

//Greatest Common Divisor
//递归
long long GCD(const long long &a,const long long &b) {
    return b?GCD(b,a%b):a;
}

//非递归~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
long long GCD(long long a,long long b) {
    while(b) {
        long long t=a%b;
        a=b;
        b=t;
    }
    return a;
}

//位运算~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
long long GCD(long long a,long long b) {
    long long d=1;
    while (a&&b)
        if (~a&1)
            if (~b&1) d<<=1,a>>=1,b>>=1;
            else a>>=1;
        else if (~b&1) b>>=1;
        else {
            if (a<b) a^=b,b^=a,a^=b;
            long long t=b;
            b=(a-b)>>1;
            a=t;
        }
    return d*(b?b: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' =d
a'(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')

这么说，你们晕了吗？没晕?真聪明。。。。

那就看代码吧！
 

//Extended GCD
long long EXGCD(const long long &a,const long long &b,long long &x,long long &y) {
    if (!b) {
        x=1;
        y=0;
        return a;
    }
    long long d=EXGCD(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

//main里面这么写
if (m%(d=EXGCD(a,b,x,y))) printf("No Solution\n");
else {
    x=x*(m/d) % (b/d);//防止x过大，使a*x爆long long
    X=a*x % b;
    printf("%lld\n",X<0?X+b:X);//之所以不用上面所讲的计算方法来求这个东西，是因为那么写太多模运算了。。看着烦。。而且感觉会比较慢。。
}

/**********************中国剩余定理***********************/

听到这名字觉得好光荣呃～～～
这也就是经常听到的孙子定理啦，韩信点兵啦什么的。。其实是一个东西。

问题是，求满足以下模线性方程组的 X（通常求的是X的最小正整数解）
X mod m1=r1
X mod m2=r2
...
...
...
X mod mn=rn
其中，mi之间两两互质

这么解：

首先 令 M=m1*m2*...*mn
对于每一个i,解 (M/mi)x+(mi)y=ri的解(x,y)
因为，mi间两两互质，那么M/mi和mi肯定是互质的，也就是说GCD(M/mi,mi)=1，方程必然有解。
这样，得到一个对应的Xi=(M/mi)*x
这个Xi有这么几个性质：

1、Xi mod mi = ri 
方程变形一下就知道了。
(M/mi)x+(mi)y=ri
<==>(M/mi)x=ri-(mi)y1
<==>(M/mi)x mod mi = ri
<==>Xi mod mi = ri
显然吧？真聪明。。。。（真的很显然。。。）

2、Xi mod mj = 0 (i!=j)
因为M/mi=m1*m2*..*m(i-1)*m(i+1)*...*mn，
也就是说，M/mi是由除了mi以外的m连乘所得，必然有
mi | M/mi
所以 mi | (M/mi)*x 
即    mi | Xi

于是，我们令SumX=X1+X2+...+Xn，那么这个 SumX 必然是满足条件的一个解。。稍微想想。能理解吧？
SumX的连加项中，只有一项Xi满足 Xi mod mi = ri,其他项Xj满足 mi | Xj (i!=j)。。懂了吧？真聪明。。。
因为每一个Xi都是存在的（刚才证明过了，GCD(M/mi,mi)=1)，所以SumX必然是存在的，也就是说，！！！这种mi两两互质的情况下肯定是有解的。！！！（注意哈！）

那么就构造出一个解SumX了。。
接下来，显然 X=SumX+k*M （k为整数） 也是方程组的解。（因为mi | M）
所以类似的，我们要得到X的最小正整数解，就可以像讲EXGCD一样的方法求了。
X%=M;
if (X<0) X+=M;

因为没有碰到中国剩余定理的裸题。。所以就没有程序了。。。。好吧？真聪明。。。。

不知道你注意到没有。。用中国剩余定理有一个前提条件，就是mi之间是两两互质的！知道我什么意思吧？真聪明。。。
如果mi并不满足两两互质的话。。怎么解呢？

说实话我也不知道为什么没有两两互质的情况下，除了无解的情况，中国剩余定理在解方程的时候会出现问题。。虽然我试过，当mi之间成倍数关系时，中国剩余定理没办法解，否则好像是可以解的。而且判无解，中国剩余定理是能判的。也就是 GCD(a,b) | m 不满足的时候，无解。

Next...


/**********************一般模线性方程组***********************/

同样是求这个东西。。
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
熟悉吧？

此时，因为GCD(a,b)=1不一定成立，GCD(a,b) | m 也就不一定成立。所以应该先判 若 GCD(a,b) | m 不成立，则！！！方程无解！！！。
否则，继续往下。

解出(x,y)，将k1=x反代回（*），得到X。
于是X就是这两个方程的一个特解，通解就是 X'=X+k*LCM(m1,m2)
这个式子再一变形，得 X' mod LCM(m1,m2)=X
这个方程一出来，说明我们实现了(1)(2)两个方程的合并。
令 M=LCM(m1,m2)，R=r2-r1
就可将合并后的方程记为 X mod M = R。

然后，扩展到n个方程。
用合并后的方程再来和其他的方程按这样的方式进行合并，最后就能只剩下一个方程 X mod M=R，其中 M=LCM(m1,m2,...,mn)。
那么，X便是原模线性方程组的一个特解，通解为 X'=X+k*M。

如果，要得到X的最小正整数解，就还是原来那个方法：

X%=M;
if (X<0) X+=M;

这么一来～～大功告成～～
中国剩余定理的题貌似都可以直接用这种方法来求了。。所以不用再特地敲中国剩余定理了。。。

代码～～～～～

long long EXGCD(const long long &a,const long long &b,long long &x,long long &y) {
    if (!b) {
        x=1;
        y=0;
        return a;
    }
    long long d=EXGCD(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
//Modular Linear Equations
//虽然它不是中国剩余定理。。。但是还是喜欢把它命名为CRT——Chinese Remainder Theorem
long long CRT(long long *m,long long *r,const int &n) {
    long long M=m[1],R=r[1],x,y,d;
    for (int i=2; i<=n; i++) {
        if ( (r[i]-R)%(d=EXGCD(M,m[i],x,y)) ) return -1;
        x=(r[i]-R)/d*x%(m[i]/d);//处理x过大。。
        R+=x*M;
        M=M*m[i]/d;
        R%=M;
    }
    return R>0?R:R+M;
}