使用扩展欧几里得算法解决POJ1061

前段时间我在刷题的时候碰到这样一道题目，题目是这样的：

http://poj.org/problem?id=1061

Description

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

 

我花费了不少时间去编写代码，但是提交后超时了。经过反思，我觉得应该是自己姿势水平不够。遂去学习了一下，认识到这道题目的思路应该是求解一个带模乘法的逆，方法就是扩展欧几里得算法。

 

Part Ⅰ：

一开始我是这么思考问题的，这两只青蛙跳了n次后见面有点像一个追及问题。他们有一个距离差D（distance）和速度差（每次能跳的距离不同）V（deltaV），如果青蛙A跳得步子更大那么就是A在追B，反之就是B在追A。

有了速度差，我们就可以不妨这样想。让前面一只静止且在D位置的地方，后面那只在0位置以速度V追它。那么问题就变成要求出后面那只青蛙绕着一个周长为L圆（也就是纬度线）跳多少次才能达到D这个地方。如果这时V=0那么肯定就不可能碰面了。设这个蛤蟆跳了n次，那么他跳了第n次的位置应该在

image.png这个地方。

那么我们就可以列出方程

image.png

可以转化为如下方程，t为一个整数

image.png

我们要求这个方程的整数解(n,t)。

 

Part Ⅱ：

为了求解这个方程，我先要介绍一个定理。

定理 1 如果任意整数a和b不都为0，则gcd(a,b)是a与b的线性组合集{ax+by: x,y∈Z}中的最小正元素。

gcd(a,b)为整数a,b的最大公约数。证明见注释。

也就是说方程image.png一定有整数解。

如果我们找到了这个方程的一个解image.png，就可以知道这个方程的基础解系。同时我们也知道了如果D不能被gcd(V,L)整除那么方程image.png就没有整数解了，也就Impossible了。因为根据定理1我们知道gcd(V,L)是V和L所能组合出的最小正整数。

那么问题来了如何找出这个方程的整数解呢？这就要用到扩展欧几里得算法了。

 

Part Ⅲ：

我们都知道欧几里得算法，也就是所谓的辗转相除法。它可能几种最古老的算法之一，它可以很快速求得两个整数的最大公约数。而扩展欧几里得算法可以收集欧几里得算法计算过程中的额外信息，使得我们能够找到image.png。

先上扩展欧几里得算法的伪代码

1
2
3
4
5
6
1 EXTENDED-EUCLID(a,b)
2     if b==0
3         return (a,1,0)
4     else(d',x',y')=EXTENDED-EUCLID(b,a mod b)
5         (d,x,y)=(d',y',x'-a/b*y')
6         return (d,x,y)

这段代码返回一个三元组，第一个元素是最大公约数，x,y 是满足ax+by=d的解。那么为什么x,y就是这个方程的整数解，它又是怎么得到的呢

 

 

下图演示了用EXTENDED-EUCLID计算gcd(99,78)的过程。

image.png

仔细观察这张图，从最后一行看。这个递归函数于b=0结束递归，返回了(3,1,0)，很明显这时的3就是最大公约数。为什么是1,0呢，很明显嘛这时a=3,b=0,x=1,y=0,很自然的就满足image.png这个方程，你y自然也可以取别的数也是满足方程的。

注意看第四行代码三元组保存了函数返回的(d',x',y')，如果这是来组最后一层递归的话那么返回的就是(3,1,0)，然后我们要求出这一层的(d,x,y)，使之满足image.png

通过刚才返回的(d',x',y')很显然满足image.png

然后我们知道image.png代换之

得到image.png

简化得image.png，第五行代码就是这么来的。

这样一层层递归函数逐渐收束，最终我们得到了image.png

 

Part Ⅳ：

借助扩展欧几里得算法，我们已经可以计算出image.png的一个整数解image.png。我们的目标是找到最小的正整数n，n存在的条件是image.png，原因见part Ⅱ。

image.png方程两边同乘image.png

得到

image.png

接下来就是通过这个方程的基础解系找到最小的n即可。

在最后一步我又踩了一次坑，那就是

image.png的通解不是

image.png而是image.png至于为什么就留给你去思考了。

 

下面附上拙劣的C代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#include<stdio.h>
struct triple
{
    long long GCD;
    long long x;
    long long y;
};
struct triple ex_gcd(long long a, long long b)
{
    struct triple temp;
    if (b == 0)
    {
        temp.GCD = a;
        temp.x = 1;
        temp.y = 1;
        return temp;
    }
    else
    {
        long long fuck;
        temp = ex_gcd(b, a%b);
        fuck = temp.y;
        temp.y = temp.x - a / b*temp.y;
        temp.x = fuck;
    }
    return temp;
}
int main()
{
    long long x, y, m, n, L;
    long long distance, deltav;
    scanf("%lld%lld%lld%lld%lld", &x, &y, &m, &n, &L);
    x = x > 0 ? x : L + x;
    y = y > 0 ? y : L + y;
    deltav = m > n ? m - n : n - m;
    distance = m > n ? (y > x ? y - x : L + y - x) : (x > y ? x - y : L - y + x);
    deltav %= L;//now deltav < L
    if (deltav == 0)
    {
        printf("Impossible");
    }
    else
    {
        struct triple key;
        key = ex_gcd(deltav, L);
        if (distance%key.GCD == 0)
        {
            if (key.x < 0)
                key.x += ((-key.x) / L + 2)*L;
            key.x *= (distance / key.GCD);
            key.x %= L / key.GCD;
            printf("%lld", key.x);
        }
        else
        {
            printf("Impossible");
        }
    }
}

注释：

见《算法导论》p545定理31.2

posted on 2017-07-31 20:29  xeoncdy  阅读(908)  评论(0)  编辑  收藏  举报