基础知识
矩阵乘法
一张图说明足矣:
代码实现(C)
const int N = 100; int c[N][N]; void multi(int a[][N], int b[][N], int n) { memset(c, 0, sizeof c); for(int i = 1; i <= n; i++) { for(int j = 1; j <= n; j++) { for(int k = 1; k <= n; k++) c[i][j] += a[i][k] * b[k][j]; } } return; }
另一种实现,将二层循环与三层循环对调,从而可以在第二层for循环时先判断 if(a[i][k] == 0) continue; ,对于矩阵有较多0的有一定效果。
int c[N][N]; void multi(int a[][N], int b[][N], int n) { memset(c, 0, sizeof c); for(int i = 1;i <= n; i++) { for(int k = 1; k <= n; k++) { if(a[i][k] == 0) continue; for(int j = 1; j <= n; j++) { c[i][j] += a[i][k] * b[k][j]; } } } return; }
显然矩阵乘法的复杂度是$O(n^3)$。
矩阵乘法需要A的行数与B的列数相等的(这是A*B的前提条件,可见矩阵乘法是不满足交换律的)。
矩阵快速幂只会用到方阵(除非题目是裸的矩阵乘法),矩阵快速幂都是方阵也就避免了检验相乘的前提条件,可以放心用。
快速幂
经典例题:数A的幂
我们利用分治的算法思想可以考虑如下求解一个数A的幂:
例如:X的19次方。
19的二进制为1 0 0 1 1。由$X^m * X^n = X^{m+n}$,则$X^{19} = X^{16} * X^2 * X^1$。
代码实现(C++)
double QuickPow(double x, int N) { if (N == 0) return 1; if (x == 0) return 0; if (N < 0) { x = 1 / x; N = -N; } double ans = 1.0; while(N) { if(N & 1) { ans *= x; } x *= x; N >>= 1; } return ans; }
过程解析
对于$X^19$来说,初始: ans = 1; res = x; 则
10011最后一位是1,所以是奇数:
ans = res * ans = x; res = res * res = x^2;
然后右移一位,1 0 0 1,则1001最后一位是1,所以是奇数:
ans = res * ans = x * (x^2) = x^3; res = res * res = x^2 * x^2 = x^4;
然后右移一位,1 0 0,则最后一位是0,所以当前的数为偶数:
res = res * res = x^4 * x^4 = x^8;
然后右移一位,1 0,最后一位是0,当前数是偶数:
res = res * res = x^8 * x^8 = x^16;
然后右移一位,1,最后一位是1,当前数是奇数:
ans = ans * res = (x^3) *(x^16) = x^19; res = res * res = x^32;
可以看出res = X^m,m始终是与二进制位置上的权值是相对应的。当二进制位为0时,我们只让res幂指数为2,对应下一个二进制位的权值;当二进制位为1时, ans = ans * res; ,则乘上了该乘的X幂次。
矩阵快速幂
上面介绍了整型快速幂的实现思路,而矩阵乘法可以用来优化很多线性递推式,那么将其中的乘法换成矩阵乘法,就变成了矩阵快速幂。
代码实现(C)
struct Matrix { //定义一个矩阵类型的结构体 int m[maxN][maxN]; } ans, res; //实现矩阵乘法A * B, n为方阵的阶 Matrix Mul(Matrix A, Matrix B, int n) { Matrix tmp; memset(tmp.m, 0, sizeof(tmp.m)); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { for (int k = 0; k < n; k++) tmp.m[i][j] = A.m[i][k] + B.m[k][j]; } } return tmp; } //快速幂算法,求矩阵res的N次幂 void QuickPower(int N, int n) { //整数快速幂ans初始化为1,对于矩阵乘法来说,则ans初始化为单位矩阵 //有A * E = A memset(ans.m, 0, sizeof(ans.m)); for (int i = 0; i < n; i++) ans.m[i][i] = 1; while (N) { if (N & 1) ans = Mul(ans, res); res = Mul(res, res); N = N >> 1; } }
经典例题:斐波拉契数列
众所周知:斐波那契数列递推公式为:$F[n] = F[n-1] + F[n-2]$, 由$F[0] = 0$,$F[1] = 1$,可以递推后面的所有数。
使用for循环是最简单直接的解法,然而对于像POJ 3070这样的题目,求斐波那契数列,其n可高达10亿,直接递推有其局限性:
(1)当n过大,递推次数过多,而测试时间有限时,可能会导致超时。
(2)想要打表实现随机访问根本不可能,先把斐波那契数列求到10亿,然后想去进行随机访问。题目未给出那么多内存,数组也开不到10亿。
因此另一个更快的解法是矩阵快速幂。
建立矩阵递推式,找到转移矩阵
显然,A就是转移矩阵。
构建矩阵递推的思路
$$A * F(n-1) = F(n)$$
一般$F(n)$与$F(n-1)$都是按照原始递推式来构建的,当然可以先猜一个$F(n)$,主要是利用矩阵乘法凑出矩阵$A$,第一行一般就是递推式,后面的行就是不需要的项就让与其相乘的系数为0。矩阵$A$就叫做转移矩阵(一定是常数矩阵),它能把$F(n-1)$转移到$F(n)$,然后这就是个等比数列,直接写出通项:$F(n) = A^{n-1} * F(1)$,此处$F1$叫初始矩阵。所以用一下矩阵快速幂然后乘上初始矩阵就能得到$F(n)$。例如,这里的$F(n)$就两个元素(两个位置),根据自己设置的$F(n)$对应位置就是对应的值。
示例
1. $f(n) = a * f(n-1) + b * f(n-2) + c$(a, b, c是常数)
$$left( egin{array}{ccc} a & b & 1 \ 1 & 0 & 0 \ 1 & 1 & 1\ end{array} ight) * left( egin{array}{c} f_{n-1} \ f_{n-2} \ c \ end{array} ight) = left( egin{array}{c} f_n \ f_{n-1} \ c \ end{array} ight)$$
2. $f(n) = c^n - f(n-1) $(c是常数)
$$left( egin{array}{cc} -1 & c \ 0 & c \ end{array} ight) * left( egin{array}{c} f_{n-1} \ c^{n-1} \ end{array} ight) = left( egin{array}{c} f_n \ c^n \ end{array} ight)$$
(整理自网络)
参考资料:
https://www.cnblogs.com/cmmdc/p/6936196.html
https://blog.csdn.net/wust_zzwh/article/details/52058209
class Solution {public: double myPow(double x, int n) { if (n == 0) return 1; if(x==0) return 0; long N = n; if (n < 0) { x = 1/x; N = -N; } double res = 1; //快速幂 while (N > 0) { if ((N & 1) == 1) { res = res * x; } x *= x; N >>= 1; } return res; }};