矩阵快速幂

关于什么是快速幂：快速幂代码
现在回想下问题斐波那契数列：$f(1)=1,f(2)=1$。在$n>2$时， $f(n)=f(n-1)+f(n-2)$。求$f(n)$。
朴素动态规划解法：

int fib(int n) {
    int dp[n];
    dp[0] = 1, dp[1] = 1;
    for (int i = 2; i < n; i++) {
        dp[i] = dp[i - 1] + dp[i - 2];
    }
    return dp[n - 1];
}

现在考虑$F_n$是由$F_{n-1}$和$F_{n-2}$线性变换得出。则有：
$$[F_n{F}_{n-1}]=[F_{n-1}{F}_{n-2}]\times \left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$
设矩阵$\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$为$P$。
则有
$$[F_{n+1}{F}_n]=[F_n{F}_{n-1}]\times \left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]=[F_{n-1}{F}_{n-2}]\times {\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^2=[F_{n-1}{F}_{n-2}]\times P^2$$
因为矩阵乘法满足结合律，所以我们可以直接用矩阵$P$的幂来得到$F_n$的值。而求矩阵$P^n$的时间复杂度为$O(log(n))$。
C++代码如下

const int maxn = 105;
struct Matrix{
    int n, m;
    int v[maxn][maxn];
    Matrix(int n, int m) : n(n), m(m) {}

    void init() {
        memset(v, 0, sizeof(v));  
    }

    Matrix operator* (const Matrix B) const {
        Matrix ans(n, B.m); // for ans
        ans.init();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < B.m; j++) {
                for (int k = 0; k < m; k++) {
                    ans.v[i][j] = ans.v[i][j] + v[i][k] * B.v[k][j];
                }
            }
        }
        return ans;
    }

    void print() {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                cout << v[i][j] << " ";
            }
            cout << endl;
        }
    }
};

Matrix q_pow (Matrix& A, int b) {
    Matrix ret(A.n, A.m);

    ret.init();
    for (int i = 0; i < ret.n; i++) { // 初始化E
        ret.v[i][i] = 1;
    }

    while (b) {
        if (b & 1) {
            ret = ret * A;
        }
        A = A * A;
        b >>= 1;
    }
    return ret;
}

完整代码如下:

#include <iostream>  
#include <cstdio>
#include <cstring>
using namespace std;
using ll = long long;
const int maxn = 105;
struct Matrix{
    int n, m;
    int v[maxn][maxn];
    Matrix(int n, int m) : n(n), m(m) {}

    void init() {
        memset(v, 0, sizeof(v));  
    }

    Matrix operator* (const Matrix B) const {
        Matrix ans(n, B.m); // for ans
        ans.init();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < B.m; j++) {
                for (int k = 0; k < m; k++) {
                    ans.v[i][j] = ans.v[i][j] + v[i][k] * B.v[k][j];
                }
            }
        }
        return ans;
    }

    void print() {
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                cout << v[i][j] << " ";
            }
            cout << endl;
        }
    }
};

Matrix q_pow (Matrix& A, int b) {
    Matrix ret(A.n, A.m);

    ret.init();
    for (int i = 0; i < ret.n; i++) { // 初始化E
        ret.v[i][i] = 1;
    }

    while (b) {
        if (b & 1) {
            ret = ret * A;
        }
        A = A * A;
        b >>= 1;
    }
    return ret;
}

int fib(int n) {
    int dp[n];
    dp[0] = 1, dp[1] = 1;
    for (int i = 2; i < n; i++) {
        dp[i] = dp[i - 1] + dp[i - 2];
    }
    return dp[n - 1];
}

int main()
{
    int n = 44;
    clock_t startTime, endTime;
    startTime = clock(); // cpu clock time
    cout << fib(n) << endl;
    endTime = clock();
    cout << (endTime - startTime) << endl;
    // [F(n)  F(n - 1)] = [F(n - 1)  F(n - 2)] * [1 1]
    //                                           [1 0]
    Matrix start = Matrix(1, 2);
    start.v[0][0] = 1;
    start.v[0][1] = 1;

    Matrix P = Matrix(2, 2);
    P.v[0][0] = 1;
    P.v[0][1] = 1;
    P.v[1][0] = 1;
    P.v[1][1] = 0;

    startTime = clock(); // cpu clock time
    Matrix ret = start * q_pow(P, n - 1);
    endTime = clock();
    cout << (endTime - startTime) << endl;
    ret.print();
}

这里暂且不考虑溢出问题。当$n$足够大的时候，可以看到快速幂需要的时钟周期明显缩短。