『矩阵入门和矩阵乘法加速递推』

---

矩阵

定义

矩阵是一个按照长方形排列的元素集合。简单地说，矩阵可以理解为一个二维数组，其中每一个位置都存放了一个元素。

例如，以下是一个大小为\(n*m\)的 矩阵。

\[ A= \left[ \begin{matrix} a_{11} & a_{12} & \cdots & a_{1m} \\ a_{21} & a_{22} & \cdots & a_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nm} \\ \end{matrix} \right] \]

矩阵运算

接下来，我们将定义矩阵的基本运算。

矩阵加法

假设有两个\(n*m\)的矩阵\(A,B\)，则矩阵\(C=A+B\)满足：
\[\forall\ i\in[1,n],j\in[1,m]\ \ \ \ \ C_{ij}=A_{ij}+B_{ij}\]
\(C\)也是一个大小为\(n*m\)的矩阵。

矩阵减法

假设有两个\(n*m\)的矩阵\(A,B\)，则矩阵\(C=A-B\)满足：
\[\forall\ i\in[1,n],j\in[1,m]\ \ \ \ \ C_{ij}=A_{ij}-B_{ij}\]
\(C\)也是一个大小为\(n*m\)的矩阵。

矩阵乘法

假设有一个\(n*m\)的矩阵\(A\)，一个\(m*p\)的矩阵\(B\)，则矩阵\(C=A*B\)满足：
\[\forall\ i\in[1,n],j\in[1,p]\ \ \ \ \ C_{ij}=\sum_{k=1}^{m}A_{ik}*B_{kj}\]

\(C\)是一个大小为\(n*p\)的矩阵。

矩阵乘法可能比较复杂，形象的理解，矩阵\(C\)第\(i\)行第\(j\)列的数，是由矩阵\(A\)第\(i\)行的\(m\)个数和矩阵\(B\)第\(j\)列的\(m\)个数先相乘再相加得到的。

例如
\[ \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} * \begin{bmatrix} b_{11} & b_{12} \\ b_{21} & b_{22} \end{bmatrix} = \begin{bmatrix} a_{11}b_{11}+a_{12}b_{21} & a_{11}b_{12}+a_{12}b_{22} \\ a_{21}b_{11}+a_{22}b_{21} & a_{21}b_{12}+a_{22}b_{22} \end{bmatrix} \]
值得注意的是，矩阵乘法满足结合律，分配律，但不满足交换律。即\(A*B=B*A\)，\(A*(B+C)=A*B+A*C\)，但是\((A*B)*C\)不一定等于\(A*(B*C)\)。

Fibonacci

Description

In the Fibonacci integer sequence, F0 = 0, F1 = 1, and Fn = Fn − 1 + Fn − 2 for n ≥ 2. For example, the first ten terms of the Fibonacci sequence are:

0, 1, 1, 2, 3, 5, 8, 13, 21, 34, …

Given an integer n, your goal is to compute the last 4 digits of Fn.

简单题意：求第n项的斐波那契数对10000求余。

Input Format

The input test file will contain multiple test cases. Each test case consists of a single line containing n (where 0 ≤ n ≤ 1,000,000,000). The end-of-file is denoted by a single line containing the number −1.

Output Format

For each test case, print the last four digits of Fn. If the last four digits of Fn are all zeros, print ‘0’; otherwise, omit any leading zeros (i.e., print Fn mod 10000).

Sample Input

0
9
999999999
1000000000
-1

Sample Output

0
34
626
6875

解析

直接递推求斐波那契数列，时间复杂度\(O(n)\)，这里显然会\(TLE\)。

我们考虑用矩阵来求解该问题。在求解\(fib_n\)时，我们只需要调用到\(fib_{n-1}\)和\(fib_{n-2}\)，那么，我们就可以设置一个矩阵\(F(n)=[fib_{n},fib_{n+1}]\)，代表当前的状态，我们称为状态矩阵。

还是利用递推的思想，我们希望能够从矩阵\(F(n-1)\)得到矩阵\(F(n)\)。观察矩阵\(F(n-1)=[fib_{n-1},fib_{n}]\)，我们发现转以后第\(2\)列的数将作为第\(1\)列的数，而新的第\(2\)列的数是原来两列数的和，所以，我们可以构造这样一个矩阵：
\[ A= \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \]
我们发现这个常数矩阵满足\(F(n-1)*A=F(n)\)，即我们可以利用该矩阵来转移状态，所以我们称之为转移矩阵。

那么我们上述式子就是递推式：\(F(n)=F(n-1)*A\)。不难进行转换：\(F(n)=F(1)*A^{n-1}\)，已知矩阵乘法是满足结合律的，那么我们就可以用经典的快速幂算法来计算\(A^{n-1}\)的值，在\(O(2^3log_2n)\)的时间内完成递推，这就是矩阵乘法加速递推。

\(Code:\)

#include<bits/stdc++.h>
using namespace std;
#define mset(name,val) memset(name,val,sizeof name)
#define mcopy(to,from) memcpy(to,from,sizeof from)
const int Mod=10000;
long long n;
inline void mul(long long f[2],long long a[2][2])
{
    long long res[2];
    mset(res,0);
    for (int j=0;j<2;j++)
        for (int k=0;k<2;k++)
            res[j] = (res[j]+f[k]*a[k][j])%Mod;
    mcopy(f,res);
    return;
}
inline void selfmul(long long a[2][2])
{
    long long res[2][2];
    mset(res,0);
    for (int i=0;i<2;i++)
        for (int j=0;j<2;j++)
            for (int k=0;k<2;k++)
                res[i][j] = (res[i][j]+a[i][k]*a[k][j])%Mod;
    mcopy(a,res);
    return;
}
inline long long solve(void)
{
    long long f[2]={0,1};
    long long a[2][2]={{0,1},{1,1}};
    for (;n;n>>=1)
    {
        if (1&n)mul(f,a);
        selfmul(a);
    }
    return f[0];
}
int main(void)
{
    freopen("fib.in","r",stdin);
    freopen("fib.out","w",stdout);
    while ( scanf("%lld",&n) && n!=-1 )
        printf("%lld\n",solve());
    return 0;
}

总结

我们初步了解了矩阵乘法优化递推这一算法，其原理在于利用矩阵来记录递推中必要的状态，再根据题意构造出恰当的转移矩阵，利用快速幂算法进行递推优化。

可想而知，不是所有递推都能用矩阵乘法优化的，该优化适用于状态矩阵长度不大，但是递推轮数很长的递推。关于转移矩阵的定义，有如下法则：如果状态矩阵的第\(i\)个数在下一个状态中会对第\(j\)个数产生影响，则将转移矩阵中第\(i\)行，第\(j\)列的数赋值为恰当的系数。

关于时间，该算法的时间复杂度为\(O(n^3logT)\)，\(n\)为状态矩阵长度，\(T\)为递推总轮数。

---

转载于:https://www.cnblogs.com/Parsnip/p/10702561.html