矩阵乘法加速递推

最新推荐文章于 2024-05-26 17:31:07 发布

W⁡angduoyu

最新推荐文章于 2024-05-26 17:31:07 发布

阅读量468

点赞数

分类专栏：数学学习总结文章标签：矩阵线性代数算法

本文链接：https://blog.csdn.net/ecjtu2020/article/details/121049748

版权

学习总结同时被 2 个专栏收录

39 篇文章 0 订阅

订阅专栏

数学

7 篇文章 0 订阅

订阅专栏

当获得一个递推式时，要求我们求第 $N$ 项，但是这个 $\le 2×10^9$ 很大，单纯的 $O (n)$ 递推无法解决。这个时候就需要用到矩阵乘法来加速递推
矩阵乘法：
设 $A$ 是 $n * m$ 矩阵， $B$ 是 $m * p$ 矩阵，则 $C = A * B$ 是 $n * p$ 矩阵，并且 $\forall i \in [1, n], \forall j \in [1, p]: \\ \quad \\ C_{i,j} = \sum_{k=1}^{m}A_{i,k}*B_{k,j}$

例题1：斐波那契

题目链接

题意：

求斐波那契数列的第 $N$ 项， $\le 2 ×10^9$

Sol：

考虑矩阵加速
构造矩阵 $①$ ：
$\begin{pmatrix} a_n&a_{n-1} &a_{n-2}\end{pmatrix}\quad$
我们要得到矩阵 $②$ ：
$\begin{pmatrix} a_n+1&a_{n} &a_{n-1}\end{pmatrix}\quad$
发现需要构造一个 $3 \times 3$ 的矩阵，根据矩阵乘法，假设这个矩阵为 $③$ ：
$\begin{pmatrix} a & d & g \\ b & e & h\\ c & f & i\end{pmatrix}\quad$
$① \times ③$ 有：
$\begin{pmatrix} a * a_n + b * a_{n-1} + c * a_{n-2} &d * a_n + e * a_{n-1} + f * a_{n-2} & g * a_n + h * a_{n-1} + i * a_{n-2}\end{pmatrix}\quad$
在根据递推式： $a_n = a_{n-1}+a_{n-2}, a_0 = 0, a_1 = a_2 = 1$
可以得出矩阵 $③$ 为：
$\begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1\\ 0 & 0 & 0\end{pmatrix}\quad$
当我们要求第 $N$ 项时，只需求出下式，然后得到的矩阵第一行第一列就是结果： $\geq 3 \\ \quad \\ \begin{pmatrix} a_2&a_1 &a_0\end{pmatrix} * \begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1\\ 0 & 0 & 0\end{pmatrix} ^ {n - 2}$

Code：

#include <bits/stdc++.h>
using namespace std;
const int mod = 10000;
struct matrix{
    int a[5][5];
}ans, base;

matrix mul1(matrix A, matrix B)
{
    matrix C;
    for(int i = 0; i < 5; ++i)
        for(int j = 0; j < 5; ++j)
            C.a[i][j] = 0;
    for(int i = 1; i <= 1; ++i)
        for(int j = 1; j <= 3; ++j)
        {
            for(int k = 1; k <= 3; ++k)
            {
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j] % mod) % mod;
                C.a[i][j] %= mod;
            }
        }
    return C;
}
matrix mul2(matrix A, matrix B)
{
    matrix C;
    for(int i = 0; i < 5; ++i)
        for(int j = 0; j < 5; ++j)
            C.a[i][j] = 0;
    for(int i = 1; i <= 3; ++i)
        for(int j = 1; j <= 3; ++j)
        {
            for(int k = 1; k <= 3; ++k)
            {
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j] % mod) % mod;
                C.a[i][j] %= mod;
            }
        }
    return C;
}

matrix qpow(int b)
{
    matrix res = ans, a = base;
    while(b){
        if(b & 1) res = mul1(res, a);
        a = mul2(a, a);
        b >>= 1;
    }
    return res;
}

int main()
{
    int n;
    while(cin >> n && n != -1)
    {
        if(n == 0) cout << 0 << endl;
        else if(n <= 2) cout << 1 << endl;
        else 
        {
            ans.a[1][1] = 1, ans.a[1][2] = 1, ans.a[1][3] = 0;
            base.a[1][1] = 1, base.a[1][2] = 1, base.a[1][3] = 0;
            base.a[2][1] = 1, base.a[2][2] = 0, base.a[2][3] = 0;
            base.a[3][1] = base.a[3][2] = base.a[3][3] = 0;
            ans = qpow(n - 2);
            cout << ans.a[1][1] << endl;
        }
    }
}

例题2：The Digits String

题目链接

题意：

求：一个长度为 $n$ 的数字字符串，他的各个数位加起来模4等于0的数量。 $\le n \le 10^9$

solution

用 $a_i$ 表示数字 $\sim 9$ 对4为 $i$ 的数量，那么有： $a_0 = 3, a_1 = 3, a_2 = 2, a_3 = 2$
然后用数组 $f_{i,j}$ 表示长度为 $i$ 的数字字符串且个数为加起来取模为 $j$ 的数量 $\le j \le3$ ，那么容易推出递推式：
$f_{i,0} = a_0 * f_{i-1,0} + a_1 * f_{i - 1,3} + a_2 * f_{i-1,2} + a_3 * f_{i-1,1} \\ f_{i,1}=a_0 * f_{i-1,1} + a_1 * f_{i - 1,0} + a_2 * f_{i-1,3} + a_3 * f_{i-1,2} \\ f_{i,2} = a_0 * f_{i-1,2} + a_1 * f_{i - 1,1} + a_2 * f_{i-1,0} + a_3 * f_{i-1,3} \\ f_{i, 3} = a_0 * f_{i-1,3} + a_1 * f_{i - 1,2} + a_2 * f_{i-1,1} + a_3 * f_{i-1,0}$
根据 $a_i$ 写出初始矩阵：
$\begin{pmatrix} 3 & 3 & 2 & 2\end{pmatrix}\quad$
根据递推式写出转移矩阵：
$\begin{pmatrix} 3 & 3 & 2 & 2 \\ 2 & 3 & 3 & 2\\ 2 & 2 & 3 & 3 \\ 3 & 2 & 2 & 3 \end{pmatrix}\quad$
那么求第 $n$ 项只需求：
$\begin{pmatrix} 3 & 3 & 2 & 2\end{pmatrix}* \begin{pmatrix} 3 & 3 & 2 & 2 \\ 2 & 3 & 3 & 2\\ 2 & 2 & 3 & 3 \\ 3 & 2 & 2 & 3 \end{pmatrix}^{n-1}$

code:

#include <bits/stdc++.h>
using namespace std;

const int mod = 2019;

struct Mat {
    int a[10][10];
};

Mat ans;

Mat mul(Mat A, Mat B) {
    Mat c;
    for(int i = 1; i <= 4; ++ i) {
        for(int j = 1; j <= 4; ++ j)
            c.a[i][j] = 0;
    }
    for(int i = 1; i <= 1; ++ i) {
        for(int j = 1; j <= 4; ++ j) {
            for(int k = 1; k <= 4; ++ k) {
                c.a[i][j] = (c.a[i][j] + A.a[i][k] * B.a[k][j]) % 2019;
            }
        }
    }
    return c;
}
Mat mul2(Mat A, Mat B) {
    Mat c;
    for(int i = 1; i <= 4; ++ i) {
        for(int j = 1; j <= 4; ++ j)
            c.a[i][j] = 0;
    }
    for(int i = 1; i <= 4; ++ i) {
        for(int j = 1; j <= 4; ++ j) {
            for(int k = 1; k <= 4; ++ k) {
                c.a[i][j] = (c.a[i][j] + A.a[i][k] * B.a[k][j]) % 2019;
            }
        }
    }
    return c;
}

Mat qpow(int b) {
    Mat base;
    base.a[1][1] = 3; base.a[1][2] = 3; base.a[1][3] = 2; base.a[1][4] = 2;
    base.a[2][1] = 2; base.a[2][2] = 3; base.a[2][3] = 3; base.a[2][4] = 2;
    base.a[3][1] = 2; base.a[3][2] = 2; base.a[3][3] = 3; base.a[3][4] = 3;
    base.a[4][1] = 3; base.a[4][2] = 2; base.a[4][3] = 2; base.a[4][4] = 3;
    while(b) {
        if(b & 1) ans = mul(ans, base);
        b >>= 1;
        base = mul2(base, base);
    }
    return ans;
}

int main()
{
    int n;
    while(cin >> n) {
        for(int i = 1; i <= 4; ++ i)
            for(int j = 1; j <= 4; ++ j)
                    ans.a[i][j] = 0;
        ans.a[1][1] = 3; ans.a[1][2] = 3; ans.a[1][3] = 2; ans.a[1][4] = 2;
        cout << qpow(n - 1).a[1][1] << endl;
    }
}