POJ1044 Date bugs

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

题目大意：

　　与众所周知的”千年虫“类似,某些计算机上存在日期记录的bug。它们的时钟有一个年份周期，每当到达最大值时，就会自动跳回到最小值。现给定一组这样的时钟，给出它们显示的年份y[i]，周期的起始年份a[i]，周期的结束年份b[i]，求可能的真实年份的最小值。真实年份应该比所有的a[i]都大。

输入：多个测试用例组成，第一行为时钟数n，接下来的n行每行为一个时钟的信息，含三个整数y[i],a[i],b[i].含义见上面的描述。n为0时表示输入结束。

输出：输出计算得的真实年份，若找不到小于10000的真实年份，输出信息："Unkown bugs detected."，用空格隔开每个用例的输出。格式见sample。

---

Sample Input

2 
1941 1900 2000 
2005 1904 2040 
2 
1998 1900 2000 
1999 1900 2000 
0 

Sample Output

Case #1: 
The actual year is 2141. 

Case #2: 
Unknown bugs detected.

---

这题与前面的一些题目相比简单太多了。记b[i]-a[i]为t[i]表示每个时钟的周期，那么真实年份可以设为y[i] + k[i] * t[i];一开始我打算先求所有周期的最小公倍数限制试探k的范围，后来发现限制了真实年份最大值才10000,很可能比最小公倍数小很多，而且公倍数还很有可能发生溢出等情况，所以还不如直接以它作为限制界限。以y[0] + k * t[0]为真实年份的试探值，假设为Y，那么对于其他的时钟都应该满足Y=y[i] + k[i] * t[i],其中k[i]为非负整数。所以只要Y - y[i] >= 0 且 Y - y[i] 可以被t[i]整除，那么这个真实年份的试探值对于第i个时钟就是可行的。所以，从0开始试探每个k，直到找到对所有时钟都可行的Y或者Y超过10000时停止即可。

//////////////////////////////////////////////////////////////////////////
//        POJ1044 Date bugs
//        Memory: 264K        Time: 0MS
//        Language: C++        Result: Accepted
//////////////////////////////////////////////////////////////////////////

#include <iostream>

using namespace std;

int n;
int y[20];
int a[20];
int t[20];

int main (void) {
    int case_id = 0;
    while (true) {
        cin >> n;
        if (n == 0) break;
        for (int i = 0; i < n; ++i) {
            cin >> y[i] >> a[i] >> t[i];
            t[i] -= a[i];    //计算出周期即可，b[i]以后不会再用到
        }
        int k = 0;
        //依次检验每个k是否可行
        while ((a[0] = y[0] + k * t[0]) < 10000) {
            bool flag = true;
            for (int i = 1; i < n; ++i) {
                int p = a[0] - y[i];
                if (p < 0 || p % t[i] != 0) {
                    flag = false;
                    break;
                }
            }
            if (flag) {
                cout << "Case #" << ++case_id <<":" << endl; 
                cout << "The actual year is " << a[0] << "." << endl << endl; 
                break;
            }     
            ++k;
        }
        if (a[0] >= 10000) {
            cout << "Case #" << ++case_id <<":" << endl; 
            cout << "Unknown bugs detected." << endl << endl; 
        }
    }
    system("pause");
    return 0;
}

---

由于一开始想要计算所有周期的最小公倍数，所以顺便回顾了一下求最小公倍数的方法（以两个数为例），在下面也记录一下吧。

小学知识告诉我们(==!)：两个数的积除以两个数的最大公约数(gcd)等于它们的最小公倍数(lcm)。

所以只要求出了gcd就很容易得到lcm了。

求gcd的经典算法是欧几里德算法，或称辗转相除法。该算法依赖与下面的性质：

记两个整数a,b的最大公约数为gcd(a, b)， 那么有定理：gcd(a, b) = gcd(b, a%b)

定理证明：设 a = k * b + r, 则 r = a % b;

假设d是a和b的一个公约数, 那么d|a(表示d可以整除a)，且d|b, 而 r = a - k * b，所以显然有d|r，故d也是a % b的约数,即 d 是 b 和 a % b 的公约数。

假设d是b和a % b的一个公约数，那么d|b, 且d|r，而 a = k * b + r, 所以有d|a, 故d也是a的约数，即d是 a 和 b 的公约数。

由以上两条知 a、b 两数的公约数与 b、 a % b两数的公约数完全一致，所以它们的最大公约数也一定相等。

所以有：gcd(a, b) = gcd(b, a%b)

辗转相除法用C++描述非常简洁：

int Gcd(int a, int b) {
    if(b == 0)
        return a;
    return Gcd(b, a % b);
}

去掉递归写为迭代形式：

int Gcd(int a, int b) {
    while(b != 0) {
        int r = b;
        b = a % b;
        a = r;
    }
    return a;
}

除了传统的欧几里德算法，还有另一种求最大公约数的算法叫Stein算法。

Stein算法主要是为了解决欧几里德算法在遇到大素数需要求模时产生的问题：

”考虑现在的硬件平台，一般整数最多也就是64位，对于这样的整数，计算两个数之间的模是很简单的。对于字长为32位的平台，计算两个不超过32位的整数的模，只需要一个指令周期，而计算64位以下的整数模，也不过几个周期而已。但是对于更大的素数，这样的计算过程就不得不由用户来设计，为了计算两个超过64位的整数的模，用户也许不得不采用类似于多位数除法手算过程中的试商法，这个过程不但复杂，而且消耗了很多CPU时间。对于现代密码算法，要求计算128位以上的素数的情况比比皆是，设计这样的程序迫切希望能够抛弃除法和取模。"

所以Stein算法中只有整数的移位和加减法，便于计算机的处理。

Stein算法基于下面的结论：

gcd(a,a) = a;

gcd(ka, kb) = k gcd(a,b);

当k与b互为质数，则gcd(ka,b)=gcd(a,b)。

算法的C++实现：

int Gcd(int a, int b) {
    if(a == 0) return b;
    if(b == 0) return a;
    if(a % 2 == 0 && b % 2 == 0) return 2 * gcd(a >> 1, b >> 1);
    else if(a % 2 == 0)  return gcd(a >> 1, b);
    else if(b % 2 == 0) return gcd(a, b >> 1);
    else return gcd(abs(a - b) / 2, min(a, b));
}

对于代码的最后一行解释如下：

这种情况下a和b都为奇数，假设a>b, 设a = b + s .那么，有gcd(a, b) = gcd(s / 2, b).

因为上面证明过gcd(a,b) = gcd(b, a%b), 这里同样令a = k * b + r,  (k >= 1)，则 s=(k - 1) * b + r.

gcd(b, s) = gcd(b, (k - 1) * b + r) = gcd(b, r); 又a、b都是奇数, 那么s一定的偶数，所以gcd(b, s) = gcd(b, s/2).

最终得到a、b都为奇数时，gcd(a, b) =  gcd(abs(a - b) / 2, min(a, b)).