整數快速冪（取模）、矩陣快速冪及其應用

摘要：

　　本文主要介紹了整數快速冪、矩陣快速冪及其應用，以題為例重點展示了使用細節。

---

 

　　我們要計算一個整數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", 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", n);
        else
            printf("%lld is normal.\n", 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\n", (st.x[0][0] - 1 + mod) % mod);
    }    
    return 0;    
}

　　關於矩陣快速冪還有其他一些重要的應用，時間有限，之后再做補充。

　　關於矩陣快速冪的介紹和應用舉例就到這里，主要運用線性代數的知識，做題的時候要找到合適的遞推式，然后利用矩陣快速冪優化。