poj 2891 扩展欧几里得迭代解同余方程组

Reference: http://www.cnblogs.com/ka200812/archive/2011/09/02/2164404.html

之前说过中国剩余定理传统解法的条件是m[i]两两互质，所以这题就不能用传统解法了= =

其实还有种方法：

 

先来看只有两个式子的方程组：

c≡b1 (mod a1)

c≡b2 (mod a2)

变形得c=a1*x+b1，c=a2*x+b2

a1*x-a2*y=b2-b1

可以用扩展欧几里得求出x和y，进而求出c

 

那么多个式子呢？可以两个两个的迭代求。

比如上面两个式子求完了，求出一个c，不妨先把它记作c1

c1满足上面两式，但对于所有的式子就不一定都满足了。

因此我们把一个新同余式加入方程组：c≡c1 (mod lcm(a1,a2))　　　　　　//lcm为最小公倍数

然后依次往下迭代就行了。最后解出的c1就是最终解。

 

如何判断无解？

对于相邻的两行a1、a2、r1、r2，若 (r2-r1) mod (gcd(a1,a2))!=0说明无解

 

#include "iostream"
using namespace std;
__int64 a[1000];
__int64 r[1000];
int n;

__int64 extend_gcd(__int64 a,__int64 b,__int64 &x,__int64 &y){
    if (b==0){
        x=1;y=0;
        return a;
    }
    else{
        __int64 r=extend_gcd(b,a%b,y,x);
        y=y-x*(a/b);
        return r;
    }
}

__int64 gcd(__int64 a,__int64 b)
{
    if (b==0) return a;
    return gcd(b,a%b);
}

__int64 lcm(__int64 a,__int64 b)
{
    __int64 tm=gcd(a,b);
    if (tm==0)  return 0;
        else return (a*b/tm);
}

int main()
{
    while (cin>>n)
    {
        for (int i=1;i<=n;i++)
            cin>>a[i]>>r[i];

        __int64 x,y,x1,c1;
        bool sol=true;
        for (int i=2;i<=n;i++)
        {
            __int64 a1=a[i-1],r1=r[i-1];
            __int64 a2=a[i],r2=r[i];
            __int64 tmp=gcd(a1,a2);
            if ((r2-r1)%tmp!=0)     sol=false;      //此处有个问题= =原来我是这样写的，WA了
            __int64 tm=extend_gcd(a1,a2,x,y);       //__int64 tm=extend_gcd(a1,-a2,x,y)
            x1=x*((r2-r1)/tm);                      //but the original equation is a1*x1-a2*x2=b2-b1
            __int64 rq=a2/tm;                       //__int64 rq=-a2/tm     
            x1=(x1%rq+rq)%rq;                       //去掉a2的负号就A了，What a fuck！！！
                                                    //个人猜想是因为gcd(a,b)和gcd(a,-b)结果是一样的，
                                                    //而且此处需要的是非负解，结果就YY对了
            c1=a1*x1+r1;
            r[i]=c1;
            a[i]=lcm(a1,a2);
           // cout<<r[i]<<" "<<a[i]<<endl;
        }
        //cout<<c1<<endl;
        __int64 ans=c1;
        if (!sol)   cout<<"-1"<<endl;
            else cout<<ans<<endl;
    }

    return 0;
}

 

 至于负号那个地方有个题hdu 1576，和本题一个情况。把负号YY掉就对了= =