模线性方程组

　　模线性方程组：x=a1(Mod m1)

                               x=a2(Mod m2)

                               ......

                               x=an(Mod mn)

　　这个问题的源自《孙子算经》，其中有这样一道算术题：“今有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二，问物几何？”按照今天的话来说：一个数除以3余2，除以5余3，除以7余2，求这个数.这样的问题，也有人称为“韩信点兵”.它形成了一类问题，也就是初等数论中解同余式.这类问题的有解条件和解的方法被称为“中国剩余定理”，这是由中国人首先提出的.首先来看一个简单的例子：

① 有一个数，除以3余2，除以4余1，问这个数除以12余几？
解：除以3余2的数有：
　　2， 5， 8， 11，14， 17， 20， 23….
它们除以12的余数是：
　　2，5，8，11，2，5，8，11，….
除以4余1的数有：
　　1， 5， 9， 13， 17， 21， 25， 29，….
它们除以12的余数是：
　　1， 5， 9， 1， 5， 9，….
　　一个数除以12的余数是唯一的.上面两行余数中，只有5是共同的，因此这个数除以12的余数是5.
　　如果我们把①的问题改变一下，不求被12除的余数，而是求这个数.很明显，满足条件的数是很多的，它是 5＋12×整数，
　　整数可以取0，1，2，…，无穷无尽.事实上，我们首先找出5后，注意到12是3与4的最小公倍数，再加上12的整数倍，就都是满足条件的数.这样就是把“除以3余2，除以4余1”两个条件合并成“除以12余5”一个条件.《孙子算经》提出的问题有三个条件，我们可以先把两个条件合并成一个.然后再与第三个条件合并，就可找到答案.
②一个数除以3余2，除以5余3，除以7余2，求符合条件的最小数.
解：　　当某数被3除余1对，即写上70（因为70是5和7的倍数，是3的倍数多1），余2时即写70×2＝140，这140仍是5和7的倍数，是3的倍数余2。某数被5除余1，即写上21（因为21是3和7的倍数、5的倍数余1），余2时，则写上21×2＝42，余3时，则写上21×3＝63。某数被7除余1时，即写上15（因为15是3和5的倍数，是7的倍数余1），余2时，则写上15×2＝30。根据题意，把70×2+21×2＋15×2计算出来结果。然后减去3、5、7的最小公倍数105，一直减到少于105为止，就得到了符合题目的数：
　　70×2＋21×3+15×2-105×2＝23
　　即此数是23。

　　通过这个例子，那么就容易退出互素的模方程的一般解法了：令M=m1*m2*...*mn，wi=M/mi，则gcd(wi,mi)=1，那么有wi*p=1(Mod mi) （其中p为未知数），则有 wi*p-mi*q=1，可以用扩展欧几里得算法求出p和q，设ei=wi*p，那么x=Σei*ai%M。

　　代码如下：

LL a[N],m[N],phi[N];
int n;

LL china()
{
    int i;
    LL M=1,w,d,y,x=0;
    for(i=0;i<n;i++)M*=m[i];
    for(i=0;i<n;i++){
        w=M/m[i];
        exgcd(m[i],w,d,d,y);
        x=(x+y*a[i])%M;
    }
    return (x+M)%M;
}

　　

　　如果不互素怎么求呢？当然也可以模仿中国剩余定理的方法，两个方程依次求解为一个方程，也就是最开始的引例：除以3余2，除以4余1，转换成12的余数是5，具体的方法如下：

　　　x=a1(Mod m1)

           x=a2(Mod m2)

     则有 m1*x+a1=m2*y+a2  =>  m1*x-m2*y=a2-a1  从这里可以看出，如果方程要有解, d|(a2-a1)必须成立，这是模线性方程有解的充分必要条件。然后可以用扩展欧几里得求出x和y了，容易得出x=x*(a2-a1)/d，M=m1/gcd(m1,m2)*m2，这样一直合并递推下去就可了。

　　代码如下：

LL a[N],m[N],phi[N];
int n;

LL Modline(int n)
{
    LL d,x,y,A,M,Mod;
    A=a[n-1],M=m[n-1];
    n--;
    // m1*x-m2*y=a2-a1
    while(n--){
        exgcd(M,m[n],d,x,y);
        if((A-a[n])%d!=0){
            return -1;
        }
        Mod=m[n]/d;
        x=(x*((a[n]-A)/d)%Mod+Mod)%Mod;
        A+=M*x;
        M=M/d*m[n];
    }
    return A;
}