[关闭]
@ArrowLLL 2017-04-06T10:56:12.000000Z 字数 3846 阅读 1836

矩阵快速幂在常系数线性递推关系中的应用

数学 算法 组合数学


主页地址 :月光森林


先引入一下,知乎上有一个问题 关于斐波拉契数列的一个低级问题 。题主询问了关于求解斐波拉契数列第n项对10007取模的结果。而这个n,可以达到甚至

解法已经在排名第一的回答中给出了,主要思路就是快速幂和矩阵乘法的结合律,亦即矩阵快速幂。具体方法这里也就不再给出。但可以依托此思想,拓展出在 的时间下计算一个递推关系的第n项。

另外要说明的一点是,这种方法仅适用于常系数线性递推关系,即递推关系

中, 以及 均为常数(可以为0)的递推关系,且 的指数均为1。

这里我们先考虑递推关系为齐次递推关系的情况,即忽略常数项b(令 ), 由线性代数的知识可以知道,在递推关系

中,如果令 , 则原递推关系可以转换为矩阵问题:

其中, A矩阵为:

可以看出, 矩阵是一个 阶方阵,且第一行由递推关系等号右边的k个元素 系数构成。从第2行到第k行为行矩阵 ——第 i 列为1、其余元素均为 0 的行矩阵,并且 。同样,这个 矩阵也可以看作是一个k阶单位矩阵去掉最后一行(即第 k 行)得到的。

接下来就要运用到矩阵乘法的结合律 来迭代上述矩阵等式:

于是可以得到等式

现在来考虑非齐次情况的递推关系,即:

。实际上读者们已经可以推出,只需要将原来的初始矩阵改变一下即可,为了方便,仍然令, 递推关系可以转化为矩阵问题:

此时A矩阵变成了 阶方阵:

这里的第一行仍然是由递推关系等号右边的k个元素 的系数以及一个构成。从第2行到第k+1行为行矩阵 ——第 i 列为1、其余元素均为 0 的行矩阵,并且 。类似地,这个 矩阵也可以看作是一个(k+1)阶单位矩阵去掉倒数第二行(即第 k 行)得到的。

通过迭代,就可以得到:

这样,便可以通过矩阵快速幂算法 快速得到递推关系 的值。做一次矩阵运算的复杂度是, 快速幂的复杂度为 。由此可以得到求解 的时间复杂度是

已知一个常系数线性递推关系

求它的第n项,n的范围是[0, 1e9]。因为数值可能过大,所以输出得数对10007取模即可。

由上面的推理我们可以得到这个常系数线性递推关系的 矩阵为

由此可以得到求解 的代码:

  1. #include <iostream>
  2. #include <cstdio>
  3. #include <cstring>
  4. using namespace std;
  5. const maxn(5), mod(10007);
  6. struct matrix{
  7. int c, r;
  8. int mat[maxn][maxn];
  9. matrix(){}
  10. matrix(int _r, int _c)
  11. :r(_r), c(_c) {memset(mat, 0, sizeof(mat));}
  12. matrix operator * (const matrix &t)const
  13. {
  14. matrix ans(r, t.c);
  15. if(c != t.r) return ans;
  16. for(int i = 0; i < r; i++)
  17. for(int j = 0; j < t.c; j++)
  18. {
  19. int &temp = ans.mat[i][j];
  20. for(int k = 0; k < c; k++)
  21. temp += mat[i][k] * t.mat[k][j];
  22. }
  23. return ans;
  24. }
  25. }A, h;
  26. void init()
  27. {
  28. int _ma[maxn][maxn] =
  29. {
  30. 2, -1, 1, 0, 0,
  31. 1, 0, 0, 0, 0,
  32. 0, 0, 1, 0, 0,
  33. 0, 0, 0, 0, 0,
  34. 0, 0, 0, 0, 0
  35. };
  36. memcpy(A.mat, _ma, sizeof(_ma));
  37. A.c = A.r = 3;
  38. h.mat[0][0] = h.mat[1][0] = 1;
  39. h.mat[2][0] = 3;
  40. h.r = 3, h.c = 1;
  41. }
  42. matrix pow(matrix t, int n)
  43. {
  44. matrix ans(3, 3);
  45. for(int i = 0; i < maxn; i++)
  46. ans.mat[i][i] = 1;
  47. while(n)
  48. {
  49. if(n & 1) ans = ans * t;
  50. t = t * t;
  51. for(int i = 0; i < t.r; i++)
  52. for(int k = 0; k < t.c; k++)
  53. t.mat[i][k] %= mod;
  54. n >>= 1;
  55. }
  56. return ans;
  57. }
  58. int main()
  59. {
  60. int n, ca = 1;
  61. while(~scanf("%d", &n) && n)
  62. {
  63. if(n == 1 || n == 2)
  64. {
  65. printf("1\n");
  66. continue;
  67. }
  68. init();
  69. matrix ans = pow(A, n - 2);
  70. ans = ans * h;
  71. printf("Case %d: %d\n",ca++, ans.mat[0][0]);
  72. }
  73. return 0;
  74. }

至于矩阵快速幂的编写方法就和上面给出的例子中类似,只需要改变matrix结构体的r和c即可,万变不离其宗。

以上です。

添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注