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

时间:2023-03-09 18:53:06
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[];
__int64 r[];
int n; __int64 extend_gcd(__int64 a,__int64 b,__int64 &x,__int64 &y){
if (b==){
x=;y=;
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==) return a;
return gcd(b,a%b);
} __int64 lcm(__int64 a,__int64 b)
{
__int64 tm=gcd(a,b);
if (tm==) return ;
else return (a*b/tm);
} int main()
{
while (cin>>n)
{
for (int i=;i<=n;i++)
cin>>a[i]>>r[i]; __int64 x,y,x1,c1;
bool sol=true;
for (int i=;i<=n;i++)
{
__int64 a1=a[i-],r1=r[i-];
__int64 a2=a[i],r2=r[i];
__int64 tmp=gcd(a1,a2);
if ((r2-r1)%tmp!=) 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 ;
}

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