整数快速幂（取模）、矩阵快速幂及其应用

摘要：

　　本文主要介绍了整数快速幂、矩阵快速幂及其应用，以题为例重点展示了使用细节。

---

　　我们要计算一个整数x的n次方，即x^n，普通的方法是连乘，这里介绍一种效率非常高的计算幂运算的算法——反复平方法。

　　首先考虑加速幂运算的方法，如果n=2^k，则可以将x^n = ((x²)²)..，即只要做k次平方运算就可以求得x^n。然后由此我们可以想到，先将n表示为2的幂次之和，即x^n = 2^k1 + 2^k2 + 2^k3... ，那么 x^n = x^{2^k1} * x^{2^k2}  * x^{2^k1} ...，只需在求x^{2^i} 的同时进行计算就好了。最终得到O(logn)的计算幂运算的算法。

　　比如计算x^22 = x^16 * x^4 * x^2，其中22的二进制数是10110，也就是需要反复平方3次。代码如下：

typedef long long ll;
ll qpow(ll x, ll n) {
    ll res = 1;
    while(n) {
        if(n&1)
            res = res * x;    //如果二进制最低位为1，则乘上x^(2^i) 
        x = x * x;            //将x平方 
        n >>= 1;            　//n/2
    }
    return res;
}

　　在实际应用中有时还需要求解x^n%mod。代码如下：

typedef long long ll;
ll qpow(ll x, ll n, ll mod) {
    ll res = 1;
    while(n) {
        if(n&1)
            res = res * x % mod;    //如果二进制最低位为1，则乘上x^(2^i) 
        x = x * x % mod;            //将x平方 
        n >>= 1;                    //n/2
    }
    return res;
}

　　看一道例题：UVA 10006 Carmichael Numbers

　　判断是否是C数，需要满足以下两个条件

　　1.不是素数.

　　2.对任意的1<x<n都有x^n和x同余模n.

　　代码如下：

#include <cstdio>
#include <cmath>
typedef long long ll;

ll qpow(ll x, ll n, ll mod) {
    ll res = 1;
    while(n) {
        if(n&1)
            res = res * x % mod;
        x = x * x % mod;
        n >>= 1; 
    }
    return res;
}
bool isprime(ll x) {
    if(x == 0 || x == 1)
        return 0;
    ll k = (ll)sqrt(x);
    for(ll i = 2; i < k; i++) {
        if(x % i == 0)
            return 0;
    }    
    return 1;
}
int main()
{
    ll n;
    while(scanf("%lld", &n) == 1 && n != 0) {
        if(isprime(n)) {
            printf("%lld is normal.
", n);
            continue;
        }
        ll i;
        for(i = 2; i < n; i++) {
            if(qpow(i, n, n) != i % n)
                break;
        }
        if(i == n) 
            printf("The number %lld is a Carmichael number.
", n);
        else
            printf("%lld is normal.
", n);
    }
    return 0;
}

　　现在要求一个矩阵A的m次幂，也就是A^m，首先应该会两个矩阵的乘法，然后知道A^m的结果一定是一个同型矩阵，最后需要理解上面的整数快速幂。剩下的就是将整数换成矩阵操作。代码如下：

const int Matrix_size = 2
struct Matrix {//定义矩阵 
    int x[Matrix_size][Matrix_size]; 
}ans, res;
/*
定义矩阵乘法，A * B， 它们的都是n阶方阵 
*/ 
Matrix Mmul(Matrix A, Matrix B, int n) {
    Matrix tmp;
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            tmp.x[i][j] = 0;
        }
    } 
    
    for(int i = 1 ; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            for(int k = 1; k <= n; k++) {
                tmp.x[i][j] += A.[i][k] * B[k][j];
            }
        }
    }
    return tmp;
}
/*
求res的N次方，n是res的阶数 
*/
void qmpow(int N, int n) {
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            res.x[i][j] = i == j ? 1 : 0;
        }
    }
    
    while(N) {
        if(N&1)
            ans = Mmul(ans, res);
        res = Mmul(res, res);
        N >>= 1;
    }
}

 　　利用上面的矩阵快速幂算法可以快速的求解一个矩阵的n次幂，那么求一个矩阵的n次幂有什么用呢？

　　1.求第n项斐波那契数

　　根据斐波那契数的定义 F₀ = 0，F₁ = 1;

　　　　　　　　　           F_n = F_{n - 1} + F_{n - 2}（n>=2）.

　　可以用矩阵表示为：

　　进一步递推得到：

　　这里需要求的是右边系数矩阵的n-2次幂，然后代入前两项即可求得f(n)，也就是第n项斐波那契数。

下面看一道例题：HDU 6198 number number number

题意

先给出了斐波那契数列的定义

 F0 = 0, F1 = 1;

 Fn = F n - 1 + F n - 2.

再给出坏数的定义，给一个整数k，如果一个整数n不能有k个任意的（可重复）的斐波那契数组成，就成为是一个坏数。现给出k，问最小的坏数是多少，答案对998244353取模。

解题思路

硬暴力的方法是不行了，因为给出的k很大。先观察前几项可得：

当k = 1时，4 = 5 - 1 = F(2 * 1 + 3) - 1;

当k = 2时，12 = 13 - 1  = F(2 * 2 + 3) - 1;

问题变成了求解第2 * k + 3项斐波那契数，又因为k很大，就需要使用矩阵快速幂解决了。

根据定义我们定义前两项是1和1，系数矩阵是1，1，1，0，求第n项需要计算系数矩阵的n-2次幂，然后乘上前两项，得到F(n)和F(n - 1)。

代码如下：

#include <cstdio>
#include <cstring>
typedef long long ll;
const int mod = 998244353;
struct Matrix {
    ll x[2][2];
};
Matrix Mmul(Matrix a, Matrix b) {
    Matrix tmp;
    memset(tmp.x, 0, sizeof(tmp.x));
    for(int i = 0; i < 2; i++) {
        for(int j = 0; j < 2; j++) {
            for(int k = 0; k < 2; k++) {
                tmp.x[i][j] = (tmp.x[i][j] + a.x[i][k] * b.x[k][j] % mod) % mod;
            }
        }
    }
    return tmp;
}
Matrix Mqpow(Matrix a, ll n) {
    Matrix tmp;
    for(int i = 0; i < 2; i++) {
        for(int j = 0; j < 2; j++) {
            tmp.x[i][j] = i == j ? 1 : 0;
        }
    }
    
    while(n) {
        if(n&1)
            tmp = Mmul(tmp, a);
        a = Mmul(a, a);
        n >>= 1;
    }
    return tmp;
}
 int main()
{
    ll k;
    while(scanf("%lld", &k) != EOF) {
        Matrix st;
        st.x[0][0] = 1; st.x[0][1] = 1;
        st.x[1][0] = 1; st.x[1][1] = 0;
        
        Matrix init;
        init.x[0][0] = 1; init.x[0][1] = 0;
        init.x[1][0] = 1; init.x[1][1] = 0;
        
        st = Mqpow(st, 2 * k + 1);
        st = Mmul(st, init);
        printf("%lld
", (st.x[0][0] - 1 + mod) % mod);
    }    
    return 0;    
}

　　关于矩阵快速幂还有其他一些重要的应用，时间有限，之后再做补充。

　　关于矩阵快速幂的介绍和应用举例就到这里，主要运用线性代数的知识，做题的时候要找到合适的递推式，然后利用矩阵快速幂优化。