POJ 1061青蛙的约会题解

POJ 1061青蛙的约会题解

网上似乎有不少此题的解法。我这个post和其他人的相比主要时想说下面几点。
1. 给出一个试图不死记硬背公式的思路;
2. 谈谈暴力解为什么看起来这么诱人，却无法通过;
3. 一些细节及对此题的评价（个人观点）

原题请参考poj: http://poj.org/problem?id=1061

我认为ACM/ICPC应该给人思维的快乐，甚至带有幽默感。这是几个世纪以来数学趣题吸引人们的那种美。
人们会废寝忘食地思考一道题目，解开的一刹那会不由自主的说: 啊--啊哈--哈哈

但是这道题目丝毫没有给我这种快感，这是一道被中国式的考试(竞赛)制度八股化的题目。
首先一看此题要求无解时给出Impossible的输出，就知道一定是要运用数论中的线性同余知识了。

真正思维的快乐在一开始的建模部分。很容易写出下面的等式:
(x + a m) mod L = (y + b n) mod L

我一开始认为两只青蛙跳跃的次数是不同的，一个为a, 另外一个为b，要求出最小的a + b。如果
题目真是是这样，那就有趣多了。可惜本题主旨是要靠死记硬背。所以实际上题意是a == b。

这样的话，上面的等式就化简为求最小的a，使得:
(x + a m) mod L = (y + a n) mod L

利用线性同余进一步化成如下形式:
a x == b (mod n)
其中 a = (n - m)
     b = (x - y)
     n = L

现在的要求就是，首先判断此线性同余方程是否有解，有解的话，给出最小正数解。
这时就必须死记硬背了。题目开始无趣了。查wikipedia得知:
[1] http://en.wikipedia.org/wiki/Linear_congruence_theorem
这个线性方程有解，当且仅当gcd(a, n)能整除b

而求gcd，当然就是递归的hello world，Euclidean算法了。可惜我脑子不好，记不得了，所以自己推导一下。
若求gcd(a, b),假设a < b, 那么一定可以写成:
b = a q + r

其中q是商，r是余数。假设d = gcd(a, b); 那么显然d能整除b和a，所以d也能整除上式的r。由于r 一定小于a
所以，我们只要求 gcd(r, a)。于是问题的规模就递归化简了。再想想边界条件，显然是 a | b。这时r为0.

所以可以轻易写出gcd的Euclidean算法了:

int gcd(int a, int b) {
    return a == 0 ? b : gcd(b % a, a);
}

于是就可以判定是否有解了。接着考虑有可能出现a或者b为负数的情形，所以需要处理下输入。正好最近读了
<<短码之美>>，索性扬弃地写出了下面的代码:

int gcd(int a, int b) {
    return a == 0 ? b : gcd(b % a, a);
}

main(a, b, x, y, m, n, L) {
    scanf("%d%d%d%d%d", &x, &y, &m, &n, &L);
    a = m < n ? n - m : m - n;
    b = m < n ? x - y : y - x;
    b = (L + b) % L; // avoid negative b
    if ( b % gcd(a, L))
        printf("Impossible\n");
    else
        ...???...
}

现在的问题是...???...的部分怎么写? 第一个思路当然是暴力解了:

if ( b % gcd(a, L))
    printf("Impossible\n");
else {
    for (x = 0, n = 0; x != b; x += a, x -= (x >= L) ? L : 0, ++n);
    printf("%d\n", n);
}

Jon Bentley在Porgramming Pearls中告诉我们，取mod是个很慢的操作，所以我就用加法和减法达到了同样的目的。

用样本测试下能通过，然后上POJ，一下子失败了，结果是Wrong Answer，结果发现这个题目的输入数据不只一行，
(虽然题目说是一行)，于是把scanf等用循环包起来。

接着还是Wrong Answer.再次读题发现，输入数据很大，是个10位数。
10位数用long能装下，但是一旦做乘法就危险了。这就是网上其他的答案中要么用long long，要么用int64的原因。
可是我的暴力解法中，没有乘法，我用加法减法化简了。所以我可以用long，改成long后上POJ，这次仍然fail了。
原因是超时Exceed Time Limit。

我自己试了试，解同余方程 11 x == 2000000000 (mod 2100000000) 的确很慢，大约用了7，8秒的样子，而题目
要求是1秒内计算出答案。

到了这里，明白这道题目就是拒绝暴力解法的。剩下的唯一思路就是快速解线性同余方程了。到这里彻底看透了考察死记硬背
的无聊用意。再次查看wikipedia，看到了Extended Euclidean algorithm:
[2] http://en.wikipedia.org/wiki/Linear_congruence_theorem#Solving_a_linear_congruence

也就是说，通过Extended Euclidean算法可以得到 a x + b y = gcd(a, b)的线性组合x和y。
这样就可以直接得到特解x0 = x b / d，其中d = gcd(a, b)。其他的解都可以表达为:
x0 + k n / d的形式，其中k为整数。

到这里已经很无趣了，实现Extended Euclidean就可以了。回想起当年看CLRS算法导论的时候，讲解RSA算法时，说过
这个扩展Euclidean。伪代码如下:

function egcd(a, b)
  if b == 0
    return (a, 1, 0)
  (d', x', y') = egcd(b, a mod b)
  return (d', y', x' - a / b * y')

翻译成C, 然后上POJ，又失败了。这次是wrong answer，这是意识到Extended Euclidean返回的x可能是负数。例如:
12 x == 20 (mod 28), 不加以处理的话，结果是-10. 于是稍作调整得到最终的程序。

翻译为C就是:

typedef long long BIG;

BIG egcd(BIG a, BIG b, BIG *x, BIG *y) {
    BIG d, x1, y1;
    if (b == 0) {
        d = a;
        *x = 1;
        *y = 0;
    }
    else {
        d = egcd(b, a % b, &x1, &y1);
        *x = y1;
        *y = x1 - a / b * y1;
    }
    return d;
}

main() {
    BIG a, b, x, y, m, n, L, d;
    while(scanf("%lld%lld%lld%lld%lld", &x, &y, &m, &n, &L) != -1) {
        a = m < n ? n - m : m - n;
        b = m < n ? x - y : y - x;
        b = (L + b) % L; // avoid negative b
        d = egcd(a, L, &x, &y);
        if ( b % d)
            printf("Impossible\n");
        else {
            for(x *= b / d, L /= d; x < 0; x += L); /*hanldle negative answer case. e.g. 12 x == 20 (mod 28)*/
            printf("%lld\n", x % L);
        }
    }
}

上POJ，终于AC了。

最后我再次认为这倒题目很无趣，整个题目的设计，尤其是细节，都强迫人走同一条道路。
由于POJ不支持Python, Scheme等(这点ZOJ就好多了)，必须使用64位整数，然后就是死记硬背扩展Euclidean算法。
好在这道题目是中文的，免得拿到世界大赛上贻笑大方了。

评论

1 楼 liuxinyu95 2013-02-28  

补充一下Extended Euclidean algorithm的推导。只要知道最基本的一条，就可以自己给出扩展Euclidean算法。

gcd(a, b)求出a, b的最大公约数，ext-gcd(a, b) 除了给出最大公约数，还给出线性组合使得：

a x + b y = gcd(a, b)

定义(d, x, y) = ext-gcd(a, b)
其中d = gcd(a, b)， x, y是上述的线性组合系数。

我们按照推导gcd的方法，通过递归减小解题规模。不是一般性，假设 a < b，于是

b = a q + r

其中q是商，r是余数。由

于d是a, b的最大公约数，故而d整除a与b。
所以上述等式中，必然有d整除r， 由于余数r是比a小的数
所以问题规模可以降低为求a与r的最大公约数
令： (d, x', y') = ext-gcd(r, a)
根据ext-gcd的定义，有：

d = x' r + y' a

从前面的 b = a q + r可得到： r = b - a q，将其带入上式，有

d = x' (b - a q) + y' a

整理得：
d = (y' - x' q) a + x' b


这正好是a 与 b的线性组合，所以有：
x = y' - x' (b / d)
y = x'

现在递归的case都推导好了。在弄好edge case就完全了。
edge case时，a = 0，此时b就是最大公约数
也就是 gcd(a, b) = b = 0 a + 1 b。

现在把上面的合并得到最终的Extended Euclidean algorithm:

ext-gcd(a, b)
  if a == 0
    return (b, 0, 1)
  (d, x', y') = ext-gcd(b mod a, a)
  return (d, y' - x' (b / d), x')

完。

后面我抽空看看能否给出线性同余方程解的推导。