矩阵快速幂学习笔记

快速幂

对于计算 xyx^y 只需要计算 xy2×xy2x^{\frac{y}{2}} \times x^{\frac{y}{2}} ，如果 yy 是奇数只需要再乘 xx ，这样就可以把 O(n)O(n) 的时间复杂度降为 O(logn)O(logn) 。

递归版：

1
2
3
4
5
6
7
8
int fastPower(int x, int y, int p) {
    if(y == 0) return 1;
    int z = fastPower(x, y >> 1, p);
    z = 1ll * z * z % p;
    if(y % 2 == 1) z = 1ll * z * x % p;
    return z;
}

迭代版：

1
2
3
4
5
6
7
8
9
10
11
int fastPower(int x,int y,int p)
{
    int ans = 1;
    while(y) {
        if(y & 1) ans = 1ll * ans * x % p;
        x = 1ll * x * x % p;
        y >>= 1;
    }
    return ans;
}

矩阵乘法

条件：第一个矩阵的列数等于第二个矩阵的行数。

得到的矩阵的行数等于第一个矩阵的行数，列数等于第二个矩阵的列数。

计算得到的矩阵第 ii 行第 jj 列的方法为取出第一个矩阵的第 ii 行与第二个矩阵的第 jj 列对应位置相乘再求和。

这里通过重载运算符简化计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
struct Matrix {
    int n, m;
    int z[110][110];
    Matrix() {
        n = m = 0;
        memset(z, 0, sizeof(z));
    }
};

Matrix operator* (const Matrix& a, const Matrix& b) {
    Matrix c;
    c.n = a.n;
    c.m = b.m;
    for(int i = 1; i <= c.n; i++)
        for(int j = 1; j <= c.m; j++) 
            for(int k = 1; k <= a.m; k++) 
            c.z[i][j] = c.z[i][j] + a.z[i][k] * b.z[k][j];
    return c;
}

矩阵快速幂

以计算斐波那契数列为例，通过递推式求第 nn 个数的时间复杂度显然是 o(n)o(n) 的，如果构造一个矩阵 MM 使得 [fifi−1]\begin{bmatrix} f_i & f_{i - 1} \end{bmatrix} 与 MM 相乘能得到 [fi+1fi]\begin{bmatrix} f_{i + 1} & f_i \end{bmatrix} ，则 fn=f1×Mn−1f_n = f_1 \times M^{n-1}，通过矩阵乘法的法则可以得到矩阵 MM 。

[fi+1fi]=[fifi−1][1110]\begin{bmatrix} f_{i + 1} & f_i \end{bmatrix} = \begin{bmatrix} f_i & f_{i - 1} \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}

所以可以通过矩阵的快速幂来优化递推。

由于已经重载了矩阵的乘号，所以只需要更改部分代码即可。

递归版：

1
2
3
4
5
6
7
8
9
10
11
12
Matrix fastPower(Matrix x, int y) {
    if(y == 0) {
        Matrix z;
        z.n = z.m = x.n;
        for(int i = 1; i <= z.n; i++) z.z[i][i] = 1;
        return z;
    }
    Matrix z = fastPower(x, y >> 1);
    z = z * z;
    if(y % 2 == 1) z = z * x;
    return z;
}

迭代版：

1
2
3
4
5
6
7
8
9
10
11
Matrix fastPower(Matrix x, long long y) {
    Matrix ans;
    ans.n = ans.m = x.n;
    for(int i = 1; i <= x.n; i++) ans.z[i][i] = 1;
    while(y) {
        if(y & 1) ans = ans * x;
        x = x * x;
        y >>= 1;
    }
    return ans;
}

矩阵乘法优化动态规划

所有能用矩阵优化的状态转移方程均为这样的形式：

f[i][j]=Σkf[i−1][k]⋅M[k][j]f[i][j] = \Sigma_kf[i - 1][k] \cdot M[k][j]

  • 一定从 f[i - 1] 转移到 f[i]。

  • 转移系数M[k][j]与i 无关。

设 ii 从 1 到 nn ，jj 与 kk 从 1 到 mm，则时间复杂度为 O(nm2)O(nm^2)。

对于上面的式子，如果将f[i]视为 1×m1 \times m 的矩阵，里面存的是 [f[i][1],f[i][2],⋯ ,f[i][m]][ f[i][1], f[i][2], \cdots, f[i][m] ]，将M视为 m×mm \times m 的矩阵，则 f[i]=f[i−1]×Mf[i] = f[i - 1] \times M ，所以 f[n]=f[0]×Mnf[n] = f[0] \times M^n ，时间复杂度变为了 O(m3logn)O(m^3logn) 。