SuperGCD(SDOI2009)

题意：

给定两个大整数 a,b $(a,b<=10^{10000})$ , 求 GCD(a,b)

分析：

高精度取模？

抱着这样的思路，我思考了一会，发现貌似行不通。

如果我们想像普通 gcd 一样取模+递归，会遇到很多麻烦：

• 首先，我们不可能实现高精/高精，这样无论代码量还是复杂度都在我们考虑的\\范围之外，如果用高精度乘法二分逼近除法结果，复杂度至少为 O(nlog^3n)，还是用 FFT 优化过的，时间复杂度和代码量都不是很理想。
• 其次，我们不能把高精度整数作为参数递归，这样绝毙爆栈，好歹也要改成” 引用”(即 C++的&)。
  所以呢？
  这让我想起了一个故事：

毕达哥拉斯学生之死

古希腊时代的伟大数学家 Pythagoras(毕达哥拉斯) 曾因为他的学生发现了根号二不可用分数表示而一怒之下将其扔进海里。追寻这个故事的缘由，原来在那个时候，人们并不会用 Euclid(欧几里得) 算法。人们在求两个数的最大公约数时，将两个数用绳子表示，每次将长的绳子减去短绳子的长度，如此反复，最终总会有一个绳子长度减为 0, 那么余下的绳子的长度就是这两个数的最大公约数。事实的检验告诉他，任何两个长度的绳子都可以反复相减使得其中一个长度变为 0, 经过推理，他断言，任何两个数都有最大公约数，因此，所有数都可以表示为另一个数的分数倍。也就是说，世界上的所有数都是可用分数表示的。

虽然说他忽略了一个事实，即绳子的精度远远不足，在相减的过程中有取整的嫌疑。但是这个故事依旧给了我们灵感：当取模失效的时候，我们还可以用减法。

然而，像毕达哥拉斯一样不停让两个数相减依旧不是什么好的办法，因为当其中一个数非常小的时候，时间复杂度将飙升至 O(n)。

所以，这个方案还需改进，使其至少达到 O(logn) 的复杂度。当然，我们还需要考虑到算法的实现，尽量用最简单的代码达到效果。

如何做到呢？我们还是需要比毕达哥拉斯聪明一点，在进行减法的同时加入一点聪明的除法。

一点点除法就够了：除以 2！

• 如果 a 和 b 同时都是偶数，我们把它们同时除以 2, 此时它们的 gcd 含有 2 这个因子，因此我们最后将答案再扩大为两倍即可。
• 直到其中一个成为奇数，或者其中一个在相减时变成了偶数，我们把它除以 2。因为此时它们的 gcd 一定不含 2 这个因子。
• 接着，我们像毕达哥拉斯一样把两数相减，继续重复这个过程。

由于每次有一个数为偶数，我们就会把它除以 2; 两个数都不是偶数的时候，相减之后又会有一个偶数。

所以这两个数在不停缩小为原来的 1/2，这样，我们就获得了一个 O(logn) 的算法。

一个看起来还不错的算法

如上的算法，只需要简单的高精度加法和高精乘 2，高精除 2 就可以实现了。而且时间复杂度非常优秀。

    #include <iostream>
    #include <cstring>
    #include <cstdio>
    #define MX 10001
    #define MB 10000
    #define MLEN 4
    using namespace std;

    int ten[6]={1,10,100,1000,10000,100000};

    typedef struct BI_t
    {
        int x[MX],len;
        BI_t operator - (const BI_t &t)
        {
            BI_t c;
            int l=max(len,t.len);
            memset(&c,0,sizeof(c));
            for(int i=1;i<=l;i++)c.x[i]=x[i]-t.x[i];
            for(int i=1;i<=l;i++)if(c.x[i]<0)c.x[i+1]--,c.x[i]+=MB;
            for(int i=l;i>=1;i--)if(c.x[i]){c.len=i;break;}
            return c;
        }
        void input()
        {
            memset(&x,0,sizeof(x));
            int slen;
            char str[MX];
            scanf("%s",str);
            slen=strlen(str);
            for(int i=0;i<slen;i++)x[i/MLEN+1]+=(str[slen-i-1]-'0')*ten[i%MLEN];
            len=(slen-1)/MLEN+1;
        }
        void out()
        {
            printf("%d",x[len]);
            for(int i=len-1;i>=1;i--)printf("%04d",x[i]);
        }
        bool operator < (const BI_t t)const
        {
            if(len!=t.len)return len<t.len;
            for(int i=len;i>=1;i--)
                if(x[i]!=t.x[i])
                    return x[i]<t.x[i];
            return 0;
        }
    }BI;

    void div2(BI &a)
    {
        if(a.x[a.len]==1)a.x[a.len]=0,a.x[a.len-1]+=MB,a.len--;
        for(int i=a.len;i>=1;i--)
        {
            a.x[i-1]+=(a.x[i]&1)*MB;
            a.x[i]/=2;
        }
    }
    void mul2(BI &a)
    {
        for(int i=1;i<=a.len;i++)a.x[i]*=2;
        for(int i=1;i<=a.len;i++)if(a.x[i]>=MB)a.x[i+1]++,a.x[i]%=MB;
        if(a.x[a.len+1])a.len++;
    }

    BI gcd(BI a,BI b)
    {
        int t=0;
        bool d1,d2;
        if(a<b)swap(a,b);
        while(b.len>1||b.x[1]!=0)
        {
            d1=!(a.x[1]&1);
            d2=!(b.x[1]&1);
            if(d1&&d2)div2(a),div2(b),t++;
            else if(d1&&!d2)div2(a);
            else if(!d1&&d2)div2(b);
            else a=a-b;
            if(a<b)swap(a,b);
        }
        while(t)mul2(a),t--;
        return a;
    }

    int main()
    {
        BI a,b;
        a.input(),b.input();
        gcd(a,b).out();
        return 0;
    }