原理

对于 $a^p$，将 $p$ 分解为 $2$ 的幂之和

• $p = \sum_{i=0}^k 2^i$
• $a^p = a^{\sum_{i=0}^k 2^i} = \prod_{i=0}^k a^{2^i}$

因此每次计算 $a^{2^i}$ 并累乘即可

步骤

求 $a^p$，迭代 $i$ 轮

• 取 $p$ 的二进制第 $i$ 位
• 如果为 $1$，则累乘 $a^{2^i}$

代码

int qpow(int a, int p){
    int res = 1;
    // p 不为0表示循坏次数小于 p 的二进制位数
    while(p){
        // p 的当前二进制位为1说明有 2^i 这一项，要乘上 a^(2^i) 这一项
        if(p & 1) res *= a;
        // 保持 a 为 a^(2^i)，每次循环 a 都要平方一次
        a *= a;
        // 右移，让 p 的下一位成为当前位，继续判断
        p >>= 1;
    }
    return res;
}

扩展

快速乘

对于 $a \times b$，将 $b$ 分解为 $2$ 的幂之和

• $b = \sum_{i=0}^k 2^i$
• $a \times b = a \times \sum_{i=0}^k 2^i = \sum_{i=0}^k a \times 2^i$

因此每次计算 $a \times 2^i$ 并累加即可

int qmul(int a, int b){
    int res = 0;
    // b 的当前二进制位为1说明有 2^i 这一项，要加上 a * 2^i 这一项
    while(b > 0){
        if(b & 1) res += a;
        // 保持 a 为 a * 2^i，每次循环 a 都要翻倍一次
        a *= 2;
        // 右移，让 b 的下一位成为当前位，继续判断
        b >>= 1;
    }
    return res;
}

矩阵快速幂

数乘换为矩阵乘法即可

• $AB = C$
• $C_{ij} = \sum_{k=1}^n A_{ik}B_{kj}$

矩阵类参考

Matrix qpow(Matrix A, int p){
    // 结果矩阵，初始化为单位矩阵
    Matrix res = Matrix(A.rows, true);
    // 要求是方阵
    if(A.rows != A.cols) return res;
    // 快速幂模板
    while(p > 0){
        // 如果当前二进制位为1，则乘上当前矩阵
        if(p & 1) res *= A;
        // 当前矩阵平方，二进制位右移
        A *= A;
        p >>= 1;
    }
    return res;
}

高精度快速幂

数乘换为高精度乘法即可

高精度参考

vector<int> qpow(vector<int> a, int p){
    vector<int> res(1, 1);
    while(p > 0){
        if(p & 1) res = mul(res, a);
        a = mul(a, a);
        p >>= 1;
    }
    return res;
}

线性递推

• 转移矩阵：将线性递推写成线性方程组，用矩阵乘法表示
• 状态向量：记录某一时刻状态的所有必要信息的向量
• 构造目标： 新状态向量 = 转移矩阵 * 旧状态向量

示例

• 状态转移方程： $a_i = a_{i-1} + a_{i-3} \space (i <= 3时a_i = 1)$
• 转移矩阵：

$$\begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}$$

• 状态向量：

$$\begin{bmatrix} a_i \\ a_{i-1} \\ a_{i-2} \\ \end{bmatrix}$$

• 构造目标：

$$\begin{pmatrix} a_i \\ a_{i-1} \\ a_{i-2} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{pmatrix} \cdot \begin{pmatrix} a_{i-1} \\ a_{i-2} \\ a_{i-3} \\ \end{pmatrix} \\ \quad\quad\quad\quad = \begin{pmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{pmatrix} ^ {i-3} \cdot \begin{pmatrix} a_{3} \\ a_{2} \\ a_{1} \\ \end{pmatrix}$$

// 转移方程
ll getK(ll k){
   if(k <= 3) return 1;
   // 构造转移矩阵和初始状态向量
   Matrix trans(3, 3), orig(3, 1, 1), res;
   trans.m[0][0] = 1;trans.m[0][1] = 0;trans.m[0][2] = 1;
   trans.m[1][0] = 1;trans.m[1][1] = 0;trans.m[1][2] = 0;
   trans.m[2][0] = 0;trans.m[2][1] = 1;trans.m[2][2] = 0;
   // 快速幂求解
   res = trans.qpow(k-3) * orig;
   // res = |a[k] a[k-1] a[k-2]｜^T
   // 取第一行第一列即为a[k]
   return res.m[0][0] % mod;
}

附录

矩阵类

struct Matrix{
    // 行数列数和矩阵数组
    int rows, cols;
    vector<vector<ll>> mat;
    // r行c列的矩阵，全员初始化(默认0)
    Matrix(int r, int c, ll init=0):
        rows(r), cols(c), mat(r, vector<ll>(c, init)) {}
    // n阶方阵，可为单位矩阵(默认非单位矩阵)
    Matrix(int n, bool identity=false): Matrix(n, n) {
        if(identity) for(int i = 0;i < n;i++) mat[i][i] = 1;
    }
    Matrix(): Matrix(2, 2){}
    // 输入输出函数
    void in(){
        for(int i = 0; i < rows; i++)
            for(int j = 0; j < cols; j++) cin >> mat[i][j];
    }
    void out(){
        for(int i = 0; i < rows; i++){
            for(int j = 0; j < cols; j++) cout << mat[i][j] << " ";
            cout << endl;
        }
    }
    // 矩阵乘法
    Matrix operator*(const Matrix b) const{
        Matrix res(rows, b.cols);
        // 要求前列数等于后行数
        if(cols != b.rows) return res;
        // 遍历前矩阵行
        for(int i = 0;i < rows;i++)
            // 遍历后矩阵列
            for(int j = 0;j < b.cols;j++)
                // 遍历前矩阵列/后矩阵行，进行点积
                for(int k = 0;k < cols;k++)
                    res.mat[i][j] += mat[i][k] * b.mat[k][j];
        return res;
    }
    Matrix& operator*=(const Matrix &b){
        return *this = *this * b;
    }
};

高精度

// 输入输出
void input(vector<int> &a){
    string s;
    cin >> s;
    for(int i = s.size() - 1;i >= 0;i--) a.push_back(s[i] - '0');
}
void output(vector<int> a){
    for(int i = a.size() - 1;i >= 0;i--) cout << a[i];
    cout << endl;
}
// 高精乘模版
vector<int> mul(vector<int> a, vector<int> b){
    int len = a.size() + b.size();
    vector<int> res(len, 0);
    for(int i = 0;i < a.size();i++){
        for(int j = 0;j < b.size();j++){
            res[i+j] += a[i] * b[j];
            res[i+j+1] += res[i+j] / 10;
            res[i+j] %= 10;
        }
    }
    while(!res.empty() && res.back() == 0) res.pop_back();
    return res;
}