数学问题（一）

1. 辗转相除法/欧几里得算法

    用辗转相除法求两个整数的最大公约数。记 gcd(a,b) 为两个数a和b的最大公约数。辗转相除法的理论依据为： gcd(a, b) = gcd(b, a % b). 
    因为设t为a和b的最大公约数，则 a = mt, b = nt, m和n互斥，a = k*b ＋ a%b， k = a/b. 那么 mt = k*nt +a%b，于是a%b可以被t整除，则t为b和a%b的公约数，设 a%b = rt. 
    t = gcd(b, a%b) <==> n 和 r 互斥 
由mt = k*nt + rt， 若 n和r不互斥，则设 n = ps, r = qs, s != 1,于是得到 m 可以被s整除，即 m和n不互斥，与之前的结论矛盾！ 
因此 gcd(a, b) = gcd(b, a%b)

算法实现

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

2. 扩展辗转相除法

    整数a和b，最大公约数为gcd(a,b)，一定有 ax + by = gcd(a, b)。如何求出x和y呢？ 
     设 bx' + a%b*y' = gcd(b, a%b)，根据欧几里得定理有bx' + a%b*y' = gcd(b, a%b) = gcd(a, b) = ax + by，则有bx' + a%b*y' = ax + by。 
由于 a%b = a - a / b，因此 bx' + (a - a/b)y' = ax + by，即 ay' + b(x' - a/b*y') = ax + by，则有 
x = y' 
y = x' - a/b*y' 
于是，可以在求得满足bx' + a%b*y' = gcd(b, a%b)的x'和y'之后再继续求出x和y，当b = 0时， x = 1， y = 0即可返回。

算法实现

//返回a和b的最大公约数，同时返回时的x和y满足 ax + by = gcd(a, b)
int extgcd(int a, int b, int& x, int& y){
    if (b == 0){
        x = 1;
        y = 0;
        return a;
    }
    int d = extgcd(b, a % b, y, x);
    y -= (a/b)*x;
    return d;
}

3. 模运算

    注意在某些环境下，若a为负数，则 a%m的结果也是负数，因此若想要得到a%m在[0, m-1)范围内的值，可以使用 (m + a % m) %m。 
    a和b对m同余，可以表示为 a = b(mod m)。假设a=c(mod m), b = d(mod m)，那么有以下基本规律： 
a + b = c + d(mod m) 
a - b = c - d(mod m) 
a x b = c x d(mod m)

4. 快速幂算法

    求一个整数x的n次方 xn，若n表示为2进制为n=ktkt−1....k1k0(b)，n=2kii1+2ki2+... 
xn=x(2ki1+2ki2+...) = x2ki1∗x2ki2∗x2ki3∗... 
于是，可以按照幂次从低到高依次计算出来x2ki，然后在n的二进制表示中若ki为1，则最后的结果需要乘以 x2ki。

long long mod_pow(long long x, long long n, long long m){
    long long res = 1;
    while(n >  0){
        if (n & 1)
            res = res * x % m;
        x = x * x % m;
        n >= 1;
    }
    return res;
}

5. 矩阵的幂

typedef vector<int> vec;
typedef vector<vec> mat;
typedef long long int ll;
const int M = 10000;

//计算 A*B
mat mul(mat& A, mat& B){
    mat C(A.size(), vec(B[0].size()));
    for(int i = 0; i < A.size(); i ++){
        for(int k = 0; k < B.size(); k ++){
            for(int j = 0; j < B[0].size(); j ++){
                C[i][j] = (C[i][j] + A[i][k]*B[k][j]) % M;
            }
        }
    }
    return C:
}
//计算A^n
mat pow(mat A, ll n){
    mat B(A.size(), vec(A.size()));
    for(int i = 0; i < A.size(); i ++){
        B[i][i] = 1;
    }
    while(n > 0){
        if (n & 1)
            B = mul(B, A);
        A = mul(A, A)
        n >= 1;
    }
    return B;
}