矩阵快速幂、矩阵加速递推

· · 个人记录

写数论文章写爽了。

矩阵快速幂

矩阵乘法

矩阵乘法要求

设 A 为 m \times p 的矩阵。B 为 p\times n 的矩阵。
矩阵乘法必须第一个矩阵的列数等于第二个矩阵的行数。

### 矩阵乘法定义 设 $A\times B=C$，则 $C_{i,j}$ 为： $$C_{i,j}=\sum^p_{k=1}A_{i,k}\times B_{k,j}$$ ### 矩阵乘法示例 假如： $$A=\begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \end{bmatrix},B=\begin{bmatrix} 1 & 2 \\ 3 & 4 \\ 5 & 6 \end{bmatrix}$$ 请你计算 $A \times B = C$ 的值，请你在草稿纸上计算出来后再翻到下面。 $\newline\textcolor{blue}{\rule{0cm}{20cm}}

我们一起计算一下：

C_{1, 2}=1\times 4+2\times 5+3\times 6=32\\ C_{2, 1}=4\times 1 + 5\times 2+6\times 3=32\\ C_{2,2}=4\times4+5\times5+6\times6=77

即：

32 & 77\end{bmatrix}

矩阵的幂运算

计算 C=\begin{bmatrix}1 & 2\\ 2 & 1\end{bmatrix}^3：

  1. 先算二次幂： \begin{bmatrix}1 & 2\\2 & 1\end{bmatrix}^2=\begin{bmatrix}1 & 2\\2 & 1\end{bmatrix}\times\begin{bmatrix}1 & 2\\2 & 1\end{bmatrix}=\begin{bmatrix}5 & 4\\4 & 5\end{bmatrix}
  2. 再算三次幂： \begin{bmatrix}5 & 4\\4 & 5\end{bmatrix}\times\begin{bmatrix}1 & 2\\2 & 1\end{bmatrix}=\begin{bmatrix}13 & 14\\14 & 13\end{bmatrix}

    单位矩阵

    单位矩阵 I=\begin{bmatrix}1 & 0\\0 & 1\end{bmatrix}，与任意矩阵相乘不改变其值。

一个 n \times n 的单位矩阵为左对角线改为 1，其他为 0 即可。

矩阵快速幂的实现

比较简单，建立一个矩阵结构体封装即可。

class matrix{
public:
    matrix()
        :n(0)
    {
        memset(c, 0, sizeof(c));
    }
    int* operator[](int t)
    {
        return c[t];
    }
    void init()//建立单位矩阵
    {
        for(int i = 1; i <= n; i++)
            c[i][i] = 1;
    }
    int size()
    {
        return n;
    }
private:
    int c[110][110];
    int n;
};
matrix operator*(matrix x, matrix y)//矩阵乘法
{
    matrix ans;
    int n = x.size();
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            for(int k = 1; k <= n; k++)
                ans[i][j] += x[i][k] * y[k][j];
}
matrix quick_pow(matrix x, int y)//矩阵快速幂
{
    matrix ans;
    ans.init();//注意这里需要是单位矩阵
    while(y)
    {
        if(y & 1) ans = ans * x;
        y >>= 1;
        x = x * x;
    }
    return ans;
}

时间复杂度为：O(n^3\times\log n )，其中 y 为指数， n 为矩阵长度。

矩阵加速递推

求斐波那契数列（P1962 斐波那契数列）：

F_{n-1}+F_{n-2} & (n > 2)\end{cases} 1\leq n \leq2^{63}

思路

构造矩阵

我们把斐波那契数列的相邻两项表示为一个矩阵 \begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}，我们希望通过 \begin{bmatrix}F_{n - 1} & F_{n - 2}\end{bmatrix} 推出 \begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}。

构造转移矩阵

我们构造一个矩阵 A，使得：

\begin{bmatrix}F_{n - 1} & F_{n - 2}\end{bmatrix}\times A=\begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}

根据递推关系：

F_{n-1}=F_{n-1}\times 1+F_{n-2}\times 0\end{cases}

观察可得：

A=\begin{bmatrix}1 & 1\\1 & 0\end{bmatrix}

递推

我们发现：

\begin{bmatrix}F_{2} & F_{1}\end{bmatrix}\times A^{n-2}=\begin{bmatrix}F_n & F_{n - 1}\end{bmatrix}

时间复杂度

矩阵快速幂为 O(n^3 \log y)，而矩阵的长度为 2，y 为题目中的 n - 2，所以时间复杂度为 O(8\log n)，就是 O(\log n)。

代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int mod = 1e9 + 7;
int n = 2;
struct matrix{
    int c[3][3];
    matrix()
    {
        memset(c, 0, sizeof(c));
    }
    void init()
    {
        for(int i = 1; i <= n; i++)
            c[i][i] = 1;
    }
};
matrix operator*(matrix x, matrix y)
{
    matrix ans;
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            for(int k = 1; k <= n; k++)
                ans.c[i][j] = (ans.c[i][j] + x.c[i][k] * y.c[k][j]) % mod;
    return ans;
}
matrix quick_pow(matrix x, int y)
{
    matrix ans;
    ans.init();
    while(y)
    {
        if(y & 1) ans = ans * x;
        y >>= 1;
        x = x * x;
    }
    return ans;
}
signed main()
{
    int m;
    cin >> m;
    if(m <= 2) cout << 1, exit(0);
    matrix a;
    a.c[1][1] = a.c[1][2] = 1;
    matrix b;
    b.c[1][1] = b.c[1][2] = b.c[2][1] = 1;
    b = quick_pow(b, m - 2);
    matrix ans = a * b;
    cout << ans.c[1][1];
    return 0;
}

作业

已知一个函数为：

f(n - 1) \times f(n - 2) \times f(n-3) & (n > 3)\end{cases}

给出 n(1\leq n\leq10^{18}) 求出 f(n)\mod 114514。

还有，欢迎来到我的博客玩哦！