通过矩阵快速幂在O(logN)时间求斐波那契数列及其推广形式

通过矩阵快速幂在O(logN)时间求斐波那契数列及其推广形式

1. 斐波那契数列

这个数列非常经典了，递推公式为：
$f(n)=f(n-1)+f(n-2)$
前几项为：
1, 2, 3, 5, 8, 13, 21, 34, 55…
如果我们要求斐波那契数列的第n项，我们一般用到的是O(N)复杂度的解法，也即构造一个长度为n的数组，第一项为1，第二项为2，接着遍历到n，按照递推公式求出每一项的值，最后返回第n项即可。因为数列每一项都可以表示为前面若干项的线性组合，所以将前面的保存下来， 避免重复计算。
实际上，斐波那契数列还有一个O(logN)复杂度的解法，通过矩阵的快速幂实现。

以下简要介绍这种思路。
因为我们知道每一项都是前面2项的线性组合，因此每推出一项就相当于对前两项组成的向量做线性变换（左乘矩阵），我们用f(k)表示数列第k项，那么我们知道f(k)=f(k-1)+f(k-2)那么
$(f(k),f(k-1){)}^T=\mathbf{M}(f(k-1),f(k-2){)}^T$
显然我们很快可以推出
$\mathbf{M}=⟮\begin{array}{r}1,1\\ 1,0\end{array}⟯$
这样根据上面的矩阵乘法，f(k)=1f(k-1)+1f(k-2),f(k-1)=1f(k-1)+0f(k-2),与公式一致。
继续推下去，则有：
$(f(k),f(k-1){)}^T={\mathbf{M}}^{k-2}(f2,f1{)}^T$
因此，求斐波那契数列的第n项，其实就是求矩阵M的n-2次方，再乘一个由第二项、第一项组成的列向量，得到的向量中的第一个。
那么矩阵的幂类似于数的幂，有O(logN)的方法。

2. 矩阵快速幂

类比于一个数的快速幂，如a的k次方，我们的思路是考虑按2进制的k分解，如：
$a^{10}=a^8\ast a^2$
也就是首先初始化结果为1，遍历k的二进制的每一位，若这一位是1，表明分解结果包含这一项，将其乘在结果上，从低位到高位（从右向左），每变化一位代表的基数就是原来的二倍，遍历时同步更新基数即可。
对应到矩阵同理：
首先我们根据矩阵乘法的定义，写一个矩阵相乘函数：

vector<vector<long long>> matMultiply(vector<vector<long long>>&A, vector<vector<long long>>&B) {//A*B
  uint64_t rA = A.size(), cA = A[0].size(), rB = B.size(), cB = B[0].size();
  if (cA != rB) {
    cout << "matrixes size error" << endl;
    return vector<vector<long long>>();
  }
  vector<vector<long long>>ans(rA, vector<long long>(cB));
  for (int r = 0; r < rA; ++r) {
    for (int c = 0; c < cB; ++c) {
      long long tmp = 0;
      for (int k = 0; k < cA; ++k) {
        tmp += A[r][k] * B[k][c];
      }
      ans[r][c] = tmp;
    }
  }
  return ans;
}

若矩阵长宽为m显然这是一个O(m^3)复杂度方法。计算过程就是矩阵乘法的定义，乘积的第i行第j列元素，等于左矩阵A第i行和右矩阵B第j列元素对应相乘之和。
接下来定义矩阵快速幂函数：

vector<vector<long long>> matFastPow(vector<vector<long long>>&A, int k) {//A^k
  if(A.empty())return vector<vector<long long>>();
  uint64_t rA = A.size(), cA = A[0].size();
  if (rA != cA) {
    cout << "matrixes size error" << endl;
    return vector<vector<long long>>();
  }
  vector<vector<long long>>ans(rA, vector<long long>(cA,0));
  for (int i = 0; i < rA; ++i) {
    ans[i][i] = 1;
  }
  for (int i = 0; i < 32; ++i) {
    if (k & 1 == 1) ans = matMultiply(ans, A);
    A = matMultiply(A, A);
    k = k >> 1;
  }
  return ans;
}

首先初始化结果为单位矩阵，接着遍历二进制k的每一位，如果是一，代表这一位对应的基数应该乘在结果上，然后更新位数与基数。
右移操作相当于移除k的最右边一位，让原来的倒数第二位变到最右边；而判断时通过k和数字1的按位与操作提取k最右边一位。
每移动一位就意味着当前基数变为原来的平方，需要调用(AA)更新。
计算过程可以认为我们使用了2logk次矩阵乘法。
因此对于斐波那契数列的第N项，这种方法的复杂度是2^3*logN,是一种O(logN)的快速解法。

3.斐波那契问题的推广

上面我们将递推公式转换成了一个系数矩阵与前几项组成向量的乘积，从而进一步推出了系数矩阵的幂乘初始值向量的形式。实际上，任何一个状态转移方程可以化简为形如斐波那契数列递推公式的dp问题，都可以采用上述的方法求解。
如果f(N)依赖他之前的m项（依赖的最前一项为f(N-m)），是这m项的线性组合，那么一定存在m*m的矩阵M，
$(f(N),f(N-1),...,f(N-m+1{)}^T={\mathbf{M}}^{N-m}(f(m),f(m-1),...,f(1){)}^T$
对于二阶的斐波那契数列形式很容易推导出M矩阵，对于更复杂的形式，可以手动推出前面若干项，自己代入计算求得M，接着利用上述矩阵快速幂方法计算即可。
对于依赖前面m项的状态转移方程，计算第N项，其时间复杂度为：
$O(m^{3}logN)$

简单举例，例如f(k)=2*f(k-1)-*f(k-3), f(1)=1,f(2)=2,f(3)=4
那么可以推出前面几项为：
1,2,4,7,12,20,33,54,…
由于依赖3项，M矩阵是3*3的，不妨设为
$\mathbf{M}=⟮\begin{array}{r}a,b,c\\ d,e,f\\ g,h,i\end{array}⟯$
则有
4a+2b+c=7
4d+2e+f=4
4g+2h+i=2
7a+4b+2c=12
7d+4e+2f=7
7g+4h+2i=4
12a+7b+4c=20
12d+7e+4f=12
12g+7h+4i=7

所以
$\mathbf{M}=⟮\begin{array}{r}2,0,-1\\ 1,0,0\\ 0,1,0\end{array}⟯$
我们输出前几项的M的幂以及乘上初始向量的结果：

可以看到得到的结果确实符合上述数列的规律。

对于状态转移方程能表示为前若干项线性组合形式的DP问题，都可以用这种方法优化到$O(m^{3}logN)$，其中m是当前项依赖的最前一项与当前一项的下标之差，换言之，f(k)能表示为f(k-1),f(k-2),…,f(k-m)的线性组合，其中系数可以为0（如上面的例子就认为f(k-2)项系数为0）。