Stein算法

        一种在高精度下求最大公约数的算法——stein算法。
        比较常规的gcd算法有辗转相除法和更相减损术，但是后者过慢，前者不利于高精度（因为有大量的除法取余运算，对高精度十分不友好），于是需要一种更高效又易于操作的算法。这便是stein算法，它是更相减损术的改进。

        算法思想如下（比如求x和y的最大公因数）：

• 如果两个数均为偶数，则两数均取半，答案为2*gcd(x/2，y/2)。
• 如果两数中有一方为偶数，一方为奇数（假如x为偶数），则偶数取半，答案为gcd(x/2，y)。
• 如果两数均为奇数，若x=y则答案即为x；否则若x>y，则答案为gcd(x-y，y)，x<y为gcd(x，y-x)。

        算法正确性很容易证明，这里略。
        重复以上过程即可得到答案。容易发现整个过程中只有除以2，乘法和减法运算，在高精度下容易操作。
        下面用洛谷P2152 SuperGCD来应用这个算法。

题目描述

        Sheng bill有着惊人的心算能力，甚至能用大脑计算出两个巨大的数的GCD（最大公约 数）！因此他经常和别人比赛计算GCD。有一天Sheng bill很嚣张地找到了你，并要求和你比赛，但是输给Sheng bill岂不是很丢脸！所以你决定写一个程序来教训他。

输入输出格式

输入格式：

        共两行： 第一行：一个数A。 第二行：一个数B。

输出格式：

        一行，表示A和B的最大公约数。

输入输出样例

Sample input

12
54

Sample output

6

说明

        对于20%的数据，0 < A , B ≤ 10¹⁸。 对于100%的数据，0 < A , B ≤ 10¹⁰⁰⁰⁰。

题解

        一道很裸的高精度gcd，注意压位（这里压8位），用stein算法就可以快速求出。
        所谓压位，就是将10进制数按10^k进制数储存在数组中（这里选k=8），从而减小时空消耗的方法。压位后的输出是一大坑点，要格外注意前导0数量的准确判定。压位后的运算过程与朴素10进位方法没有大的差别，要注意提前设置好进位base。
        这题不要用递归，否则MLE，用数组循环即可。

#include<iostream>
#include<string>

#define P 8
#define Base 100000000
using namespace std;
long long s1[10000] = {0}, s2[10000] = {0}, ans[10000] = {0}, temp[10000] = {0};

inline int change(string x) {//化字符串为数字
    int s = 0;
    for (int i = 0; i < x.length(); i++)s = s * 10 + x[i] - '0';
    return s;
}

inline bool isEven(const long long *x) {//判断偶数
    return x[1] % 2 == 0;
}

inline void div(long long *x) {//除以2
    int r = 0;
    for (int i = x[0]; i >= 1; i--) {
        if (x[i] % 2 == 0) {
            x[i] = (x[i] + r * Base) / 2;
            r = 0;
        } else x[i] = (x[i] + r * Base) / 2, r = 1;
    }
    if (x[x[0]] == 0)x[0]--;
}

inline void times(long long *x) {//x=x*2
    int r = 0;
    for (int i = 1; i <= x[0]; i++)x[i] *= 2, x[i] += r, r = x[i] / Base, x[i] %= Base;
    if (r != 0)x[++x[0]] = r;
}

inline int cmp(long long *x, long long *y) {//x<y?
    if (x[0] < y[0])return 1;
    if (x[0] > y[0])return -1;
    for (int i = x[0]; i >= 1; i--)
        if (x[i] < y[i])return 1;
        else if (x[i] > y[i])return -1;
    return 0;//=
}

inline void times2(long long *x, const long long *y) {//x=x*y
    for (int i = 1; i <= x[0]; i++) {
        for (int j = 1; j <= y[0]; j++) {
            temp[j + i - 1] += x[i] * y[j];
        }
    }
    long long r = 0;
    x[0] = y[0] + x[0] - 1;
    for (int i = 1; i <= x[0]; i++) {
        long long t = temp[i] + r;
        x[i] = t % Base;
        r = t / Base;
    }
    while (r != 0)x[++x[0]] = r % Base, r /= Base;
}

inline void sub(long long *x, long long *y) {//x=x-y
    for (int i = 1; i <= x[0]; i++) {
        if (x[i] >= y[i])x[i] -= y[i];
        else {
            x[i] = x[i] + Base - y[i];//不够借位
            x[i + 1]--;
        }
    }
    for (int i = x[0]; i >= 1; i--)//重新找右值
        if (x[i] != 0) {
            x[0] = i;
            return;
        }
}

inline void print(long long *x) {//输出x
    int te;
    for (int i = x[0]; i >= 1; i--) {
        if (i == x[0]) cout << x[i];
        else {
            te = Base / 10;
            while (te > 0)cout << x[i] / te, x[i] %= te, te /= 10;
        }
    }
    cout << endl;
}

int main() {
    string a, b;
    cin >> a >> b;
    for (int i = a.length() - 1; i >= 0; i -= P) {//压位
        if (i - P + 1 >= 0)s1[++s1[0]] = change(a.substr(i - P + 1, P));
        else s1[++s1[0]] = change(a.substr(0, i + 1));
    }
    for (int i = b.length() - 1; i >= 0; i -= P) {//压位
        if (i - P + 1 >= 0)s2[++s2[0]] = change(b.substr(i - P + 1, P));
        else s2[++s2[0]] = change(b.substr(0, i + 1));
    }
    ans[1] = ans[0] = 1;
    while (true) {//stein算法执行
        if (isEven(s1) && isEven(s2))div(s1), div(s2), times(ans);
        else if (isEven(s1))div(s1);
        else if (isEven(s2))div(s2);
        else {
            int status = cmp(s1, s2);
            if (status == 0) {
                times2(ans, s1);
                print(ans);
                return 0;
            } else if (status == 1)sub(s2, s1);
            else sub(s1, s2);
        }
    }
}

        P.S.据说Python代码就两行…Orz