poj 1061 青蛙的约会（用欧几里得解线性方程）

题目：来源于http://poj.org/problem?id=1061

转自：http://huangfeihome.blog.163.com/blog/static/208900245201277114052577/（有关欧几里得解方程的推导）

              http://hi.baidu.com/newmyl/item/dc6805253305833195f62b0a  （有关题目的推导和解法）

          两只青蛙在网上相识了，它们聊得很开心，于是觉得很有必要见一面。它们很高兴地发现它们住在同一条纬度线上，于是它们约定各自朝西跳，直到碰面为止。可是它们出发之前忘记了一件很重要的事情，既没有问清楚对方的特征，也没有约定见面的具体位置。不过青蛙们都是很乐观的，它们觉得只要一直朝着某个方向跳下去，总能碰到对方的。但是除非这两只青蛙在同一时间跳到同一点上，不然是永远都不可能碰面的。为了帮助这两只乐观的青蛙，你被要求写一个程序来判断这两只青蛙是否能够碰面，会在什么时候碰面。

我们把这两只青蛙分别叫做青蛙A和青蛙B，并且规定纬度线上东经0度处为原点，由东往西为正方向，单位长度1米，这样我们就得到了一条首尾相接的数轴。设青蛙A的出发点坐标是x，青蛙B的出发点坐标是y。青蛙A一次能跳m米，青蛙B一次能跳n米，两只青蛙跳一次所花费的时间相同。纬度线总长L米。现在要你求出它们跳了几次以后才会碰面。

Input
输入只包括一行5个整数x，y，m，n，L，其中x≠y < 2000000000，0 < m、n < 2000000000，0 < L < 2100000000。
Output
输出碰面所需要的跳跃次数，如果永远不可能碰面则输出一行"Impossible"

Sample Input

1 2 3 4 5

Sample Output

4

 
 
关于题目的推导：
 
             此题其实就是扩展欧几里德算法－求解不定方程，线性同余方程。
　             　设过s步后两青蛙相遇，则必满足以下等式：

　　　　              (x+m*s)-(y+n*s)=k*l(k=0,1,2....)

　　             稍微变一下形得：

　　　　              (n-m)*s+k*l=x-y

                    令n-m=a,k=b,x-y=c,即

　　　                　a*s+b*l=c

　　             只要上式存在整数解，则两青蛙能相遇，否则不能。

 
 
扩展欧几里德算法求出二元一次方程 a * x + b * y = gcd(a, b) 的解

     
算法用到了递归的思想。      

      设 d= gcd(a, b) = gcd(b, a % b)，结合欧几里德算法，我们首先可以把gcd(a, b) = ax + by转换为另一个方程gcd(b, a % b) = bx' + a % by'.为什么进行这种转换呢？因为通过欧几里德算法我们知道这种转换一定会停止，并且最终gcd(a,
b) = ax + by会变成gcd(a, 0) = ax + 0*y而这个方程很容易解，只要x =1，y = 0就可以满足了（当然，y可以是任何数，只是选择不同的y对应最后解集合中不同的一组解）。
      但是进行这样的转换后，很明显x' y'就不是原来的解了。根据递归的思想，假设现在我们已经求出了x'
y',剩下的关键就是如何用x' y'求出x y.我们观察gcd(b, a % b) = bx' + a % by'，因为实际上gcd(b, a % b)  =gcd(a, b)只要把右边重新写成 =ax + by就可以了，所以需要对bx'+a%by'进行变形，因为a % b = a - (a / b) * b,故bx' + a % by' = bx' + (a - a / b * b ) y' = ay' + b (x' - a / b * y')（这里的除法是程序中的向下取整的除法）。
这样便可得出：
      x = y'
      y = x' -  a / b * y'
这也就是递归回来时所用到的式子。
      这样我们就得到了ax+by=gcd(a,b)的一组解(x0, y0)，即ax0 + by0 = gcd(a,b)，实际上很容易得出一系列解
      x = x0 + b / gcd(a, b) * k
      y = y0 - a / gcd(a, b) * k  （k为任意整数）
要验证只需代入方程即可。

扩展欧几里德解一般二元一次方程a * x + b * y = c

      从上面的过程可以看到，ax + by = gcd(a ,b)一定有解，再看一般形式ax + by=c，只有当且仅当c是gcd(a, b)的整数倍时才有解，否则无解。因为如果我们两边同除gcd(a, b)，会出现这样一种情况a'x + b'y=c / gcd(a, b)，如果c不是gcd(a,b)的整数倍，那么左边是整数，右边是分数，显然无解。
      至此，到网上搜很多博客都会发现这句话“ ax+by=c的求解可以先求出ax
 + by = gcd(a, b)，然后将x y扩大c / gcd(a,b)倍就可以了“，而我自己无法证明这句话的正确性，反而觉得这句话是错的，在前面已经得出ax + by = gcd(a, b)的解是：
      x = x0 + b / gcd(a, b) * k
      y = y0 - a / gcd(a, b) * k  （k为任意整数）
如果将x，y扩大c / gcd(a, b)倍后就成了
      x = x0 * c / gcd(a, b) + b * c * k / (gcd(a, b) ^ 2)
      y = y0 * c / gcd(a, b) - a * c * k / (gcd(a, b) ^ 2)（k为任意整数）
而我解出的通解是
      x = x0 * c / gcd(a, b) + b / gcd(a, b) * k
      y = y0 * c / gcd(a, b) - a / gcd(a, b) * k  （k为任意整数）
上面两组通解明显不等价，下面的一组通解包含了上面的一组通解，通过代入原方程ax + by = c验证，下面的一组通解才是正确的。所以，在求通解的过程中，只需要在特解x0 y0上扩大c / gcd(a, b)倍就可以了。

 
代码：

#include <iostream>
using namespace std;
long long extgcd(long long a, long long b, long long &x, long long &y)
{
long long d, t;
if (b == 0) { x = 1; y = 0; return a; }
d = extgcd(b, a % b, x, y);
t = x - a / b * y; x = y; y = t; //由gcd(b, a % b) = bx' + a % by'逆推ax+by=gcd(a,b)的解
return d;
}
int main()
{
long long x, y, m, n, L, X, Y, d, r;
while (cin >> x >> y >> m >> n >> L)
{
d = extgcd(n - m, L, X, Y); r = L / d;
if ((x - y) % d) cout << "Impossible" << endl;
else cout << ((x - y) / d * X % r + r) % r << endl;
}
return 0;
}

 
 
另附本题另一篇博文的代码：

#include<iostream>
using namespace std;
long long gcd(long long a,long long b)
{
long long t;
if(b==0)
return a;
return gcd(b,a%b);
}
void exgcd(long long a,long long b,long long &m,long long  &n)
{
if(b==0)
{
m=1;
n=0;
return ;
}
exgcd(b,a%b,m,n);
long long t=m;
m=n;
n=t-a/b*n;
}
int main()
{
long long x,y,m,n,l,a,b,k1,k2,t;
while(cin>>x>>y>>m>>n>>l)
{
a=n-m;
b=l;
n=x-y;
long long r=gcd(a,b);
if(n%r)
{
cout<<"Impossible"<<endl;
continue;
}
a=a/r;
b=b/r;
n=n/r;
exgcd(a,b,k1,k2);
t=n*k1/b;//注1
k1=n*k1-t*b;
if(k1<0)
k1+=b;
cout<<k1<<endl;
}
return 0;
}

 
注1：（做废）

　 　此时方程的所有解为:x = c * k1 + b * t，令x=0可求出当x最小时的t的取值，这样求出的x可能小于0，当其小于0时加上b即可。 
 
 注1：（原注1做废）
        此时方程的所有解为：x=c*k1-b*t,x的最小的可能值是0，令x=0可求出当x最小时的t的取值，但由于x=0是可能的最小取值，实际上可能x根本取不到0，那么由计算机的取整除法可知：由 t=c*k1/b算出的t，代回x=c*k1-b*t中，求出的x可能会小于0，此时令t=t+1，求出的x必大于0；如果代回后x仍是大于等于0的，那么不需要再做修正。