矩阵快速幂

一、快速幂

实现 pow(x, n) ，即计算 x 的整数 n 次幂函数（即，xn ）。

对于上面这道题，我们最容易想到的解法就是通过循环或递归进行累乘求解

public double MyPow(double x, int n)
{
	long N = n;
	if (N == 0)
		return 1;
	double res = x;
	long num = N > 0 ? N : -N;
	for (long i = 1; i < num; i++)
	{
		res *= x;
	}

	if (N < 0)
		res = 1 / res;
	return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

然而，这种解法在面对n比较大的情况时，效率会变得十分低下。仔细想想上面的循环过程就会发现，里面进行了很多次无意义的计算。以$a^{10}$举例，通过循环一个一个往上乘，需要进行9次运算。而如果我们先计算出$a^2$（1次），然后计算$a^2\ast a^2\ast a^2\ast a^2\ast a^2$（4次），总计算量就能减少到5次。如果先计算出$a^4$（2次），再计算$a^2\ast a^4\ast a^4$（2次），总计算量就能减少到4次。

也就是说，将$a^n$拆解成$a^x\ast a^y\ast a^z...$的形式就可以降低总体运算的次数。那么如何拆解n才最高效呢？根据前面的思路，显然是每次平方所需的运算次数最少。仍然是$a^{10}$的例子，我们从$a^1$开始，每次取平方值，可以得到$[a^1,a^2,a^4,a^8]$。想要组成$a^{10}$，只需要让$a^8\ast a^2$即可。所需的计算次数也是4次。看到这里有没有觉得熟悉？我们实际上就是在把n拆成几个$2^{?}$相加，也就是类似于n的二进制转换成十进制的过程。

下面仍然以$a^{10}$为例，通过图展示一下计算过程：
我们从10的最低位开始遍历，遇到的第一位是0，对应的是$a^1$，但我们并不需要它计算结果，所以跳过。

接下来遇到的第二位是1，对应$a^2$，这是计算结果需要用到的，所以把它乘到res上。继续移动指针。

遇到的第三位是0，对应$a^4$，用不到，跳过。

遇到的第四位是1，对应$a^8$，需要用到，把它乘到res中。

res便是最终结果。将上面的思路转换成代码（未考虑特殊情况）：

private double Pow(double a, int n)
{
	double res = 1;
	while (n !=0)
	{
		// 当前位为1,计入结果
		if ((n & 1) == 1)
			res = res * a;
		// a自乘
		a *= a;
		// n右移
		n >>= 1;
	}
	return res;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

二、矩阵快速幂

矩阵快速幂的一个非常典型的应用就是斐波那契数列（原题地址）。先来看斐波那契数列的定义：
$$F(0)=0,F(1)=1$$
$$F(N)=F(N-1)+F(N-2),其中N>1.$$
我们先定义一个矩阵
[ F ( N ) F ( N − 1 ) ]

[F(N)​F(N−1)​]
现在假设有一个变换矩阵$M$，与上述矩阵存在如下关系

[ F ( N ) F ( N − 1 ) ] × M = [ F ( N + 1 ) F ( N ) ]

×M= [F(N)​F(N−1)​]×M=[F(N+1)​F(N)​]

很显然，$M$是一个$2\times 2$的矩阵，我们假设
M = [ a b c d ] M=

M=[ac​bd​]
带入上面的公式可得
[ F ( N ) F ( N − 1 ) ] × [ a b c d ] = [ a F ( N ) + c F ( N − 1 ) b F ( N ) + d F ( N − 1 ) ] = [ F ( N + 1 ) F ( N ) ]

[F(N)F(N−1)]" role="presentation" style="position: relative;">[F(N)F(N−1)]$$\left[\begin{array}{cc}F(N)& F(N-1)\end{array}\right]$$

×

[abcd]" role="presentation" style="position: relative;">[acbd]$$\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$$

=

[aF(N)+cF(N−1)bF(N)+dF(N−1)]" role="presentation" style="position: relative;">[aF(N)+cF(N−1)bF(N)+dF(N−1)]$$\left[\begin{array}{cc}aF(N)+cF(N-1)& bF(N)+dF(N-1)\end{array}\right]$$

=

[F(N+1)F(N)]" role="presentation" style="position: relative;">[F(N+1)F(N)]$$\left[\begin{array}{cc}F(N+1)& F(N)\end{array}\right]$$

[F(N)​F(N−1)​]×[ac​bd​]=[aF(N)+cF(N−1)​bF(N)+dF(N−1)​]=[F(N+1)​F(N)​]
也就是
$$aF(N)+cF(N-1)=F(N+1)$$
$$bF(N)+dF(N-1)=F(N)$$
很容易得出
M = [ 1 1 1 0 ] M =

[1110]" role="presentation" style="position: relative;">[1110]$$\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$

M=[11​10​]
即
[ F ( N ) F ( N − 1 ) ] × [ 1 1 1 0 ] = [ F ( N + 1 ) F ( N ) ]

[F(N)F(N−1)]" role="presentation" style="position: relative;">[F(N)F(N−1)]$$\left[\begin{array}{cc}F(N)& F(N-1)\end{array}\right]$$

×

[1110]" role="presentation" style="position: relative;">[1110]$$\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$

=

[F(N+1)F(N)]" role="presentation" style="position: relative;">[F(N+1)F(N)]$$\left[\begin{array}{cc}F(N+1)& F(N)\end{array}\right]$$

[F(N)​F(N−1)​]×[11​10​]=[F(N+1)​F(N)​]
根据数学归纳法可得
[ F ( N + 1 ) F ( N ) ] = [ F ( 1 ) F ( 0 ) ] × [ 1 1 1 0 ] n

[F(N+1)F(N)]" role="presentation" style="position: relative;">[F(N+1)F(N)]$$\left[\begin{array}{cc}F(N+1)& F(N)\end{array}\right]$$

=

[F(1)F(0)]" role="presentation" style="position: relative;">[F(1)F(0)]$$\left[\begin{array}{cc}F(1)& F(0)\end{array}\right]$$

×

[1110]" role="presentation" style="position: relative;">[1110]$$\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$$

^n [F(N+1)​F(N)​]=[F(1)​F(0)​]×[11​10​]n
我们只需要定义出矩阵乘法，配合前面的快速幂，即可实现在$O(logN)$时间复杂度内解出斐波那契数列。

// 自定义的矩阵类
public class MyMatrix
{
    private int[,] _matrix;
    public int RowNum { get; }
    public int ColNum { get; }

    public int this[int row, int col]
    {
        get => _matrix[row, col];
        set => _matrix[row, col] = value;
    }
    public MyMatrix(int rowNum,int colNum)
    {
        _matrix = new int[rowNum, colNum];
        RowNum = rowNum;
        ColNum = colNum;
    }
	// 矩阵乘法
    public MyMatrix Multiply(MyMatrix b)
    {
        if (ColNum != b.RowNum)
            throw new Exception("矩阵1的列数与矩阵2的行数需相同");
        MyMatrix newMatrix = new MyMatrix(RowNum, b.ColNum);
        for (int i = 0; i < RowNum; i++)
        {
            for (int j = 0; j < ColNum; j++)
            {
                for (int k = 0; k < b.ColNum; k++)
                {
                	// 根据题目要求,结果需要对1000000007取余
                    newMatrix[i, k] = (int)(((long)this[i, j] * b[j, k] + newMatrix[i, k])%1000000007);
                }
            }
        }
        return newMatrix;
    }
}
public class Solution {
	// 程序入口
    public int Fib(int n)
    {
        if (n == 0)
            return 0;
        if (n == 1)
            return 1;
        // 定义变换矩阵
        MyMatrix a = new MyMatrix(2, 2)
        {
            [0, 0] = 1,
            [0, 1] = 1,
            [1, 0] = 1,
            [1, 1] = 0
        };
        MyMatrix res = Pow(a, n);
        // 因为[f(1) f(0)] = [1,0],无需计算直接得出f(n)
        return res[0, 1];
    }

    private MyMatrix Pow(MyMatrix a, int n)
    {
        // 单位矩阵,相当于1
        MyMatrix res = new MyMatrix(2, 2)
        {
            [0, 0] = 1,
            [0, 1] = 0,
            [1, 0] = 0,
            [1, 1] = 1
        };
        while (n != 0)
        {
            if ((n & 1) == 1)
                res = res.Multiply(a);
            a = a.Multiply(a);
            n >>= 1;
        }
        return res;
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79