矩阵快速幂


矩阵快速幂


快速幂

    为了更快的进行幂运算,我们可以使用这个方法。原理太简单了,不解释。
看一下快速幂的非递归实现:

typedef long long ll;
/*
base 基数
p 指数
m 模数(看要求,决定要不要)
*/
ll fastPower(ll base, ll p, ll m)
{
    ll res = 1;
    while(p)
    {
        if(p & 1)   res = res * base % m;
        p >>= 1;
        base = base * base % m;
    }
    return res;
}

矩阵快速幂

    矩阵快速幂是为了加速矩阵的乘法运算。我在大一的线性代数课上学到了这个矩阵的相关知识,当时没觉得有什么用,现在发现...有点东西。
矩阵相乘的算法呢,百度一下有很多,我在这贴一个 OI Wiki 上的C++实现。

struct mat {
  LL a[sz][sz];
  inline mat() { memset(a, 0, sizeof a); }
  inline mat operator-(const mat& T) const {
    mat res;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) {
        res.a[i][j] = (a[i][j] - T.a[i][j]) % MOD;
      }
    return res;
  }
  inline mat operator+(const mat& T) const {
    mat res;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) {
        res.a[i][j] = (a[i][j] + T.a[i][j]) % MOD;
      }
    return res;
  }
  inline mat operator*(const mat& T) const {
    mat res;
    int r;
    for (int i = 0; i < sz; ++i)
      for (int k = 0; k < sz; ++k) {
        r = a[i][k];
        for (int j = 0; j < sz; ++j)
          res.a[i][j] += T.a[k][j] * r, res.a[i][j] %= MOD;
      }
    return res;
  }
  inline mat operator^(LL x) const {
    mat res, bas;
    for (int i = 0; i < sz; ++i) res.a[i][i] = 1;
    for (int i = 0; i < sz; ++i)
      for (int j = 0; j < sz; ++j) bas.a[i][j] = a[i][j] % MOD;
    while (x) {
      if (x & 1) res = res * bas;
      bas = bas * bas;
      x >>= 1;
    }
    return res;
  }
};

    那它有什么用呢,目前我知道的是,它可以加快一些递推公式的递推速度。
比如说呢,斐波那契数列:
f[0] = f[1] = 1,
f[n] = f[n - 1] + f[n - 2].
也就是$$a_0$$ = 1, $$a_1$$ = 1, $$a_n$$ = $$a_{n-1}$$ + $$a_{n - 2}$$。
    如果我们要求$$a_n$$的话,我们需要从$$a_0$$ = 1, $$a_1$$ = 1开始递推到$$a_n$$,时间复杂度是O(N).但是如果我们用到矩阵加速的话,时间复杂度可以小很多具体多少额也不清楚。看看下面这个式子:
$$\begin{bmatrix}{F}_{2}&{F}_{3}\end{bmatrix} = \begin{bmatrix}{F}_{1}&{F}_{2}\end{bmatrix} * \begin{bmatrix}0&1\\1&1\end{bmatrix}$$
它就是$$F_3 = F_1 + F_2$$这个过程,那就好说了。
$$\begin{bmatrix}{F}_{n}&{F}_{n + 1}\end{bmatrix} = \begin{bmatrix}{F}_{1}&{F}_{2}\end{bmatrix} * {\begin{bmatrix}0&1\\1&1\end{bmatrix}}^{n - 1}$$
现在我们要求$$F_n$$的话,只需要对右边那个矩阵进行矩阵加速,得到它的n - 1次方即可。
我写了个模板:

#define MATRIXSIZE 2
//a b矩阵相乘,结果返回给a
void matrixMultiply(int a[][MATRIXSIZE], int la, int ca, int b[][MATRIXSIZE], int lb, int cb)
{
    int i, j, k;
    int c[MATRIXSIZE][MATRIXSIZE] = {0};
    for(i = 0; i < la; i++)
        for(j = 0; j < cb; j++)
            for(k = 0; k < ca; k++)
                c[i][j] += a[i][k] * b[k][j];

    for(i = 0; i < la; i++)
        for(j = 0; j < ca; j++)
            a[i][j] = c[i][j];
}
//res中传入初始矩阵,base中传入基数矩阵,p为指数
void matrixFastPower(int base[][MATRIXSIZE], int p, int res[][MATRIXSIZE])
{
    while(p)
    {
        if(p & 1)   matrixMultiply(res, 1, MATRIXSIZE, base, MATRIXSIZE, MATRIXSIZE);   //矩阵大小手动改一下
        p >>= 1;
        matrixMultiply(base, MATRIXSIZE, MATRIXSIZE, base, MATRIXSIZE, MATRIXSIZE);
    }
}

例题

贴个题运用一下: 洛谷P1939, 点我看题
    这个题的递推公式是隔一个相加,我当时愣是没想出来怎么构建这个矩阵,实在是太菜了。看完题解之后,给了自己一个大嘴巴子。
$$\begin{bmatrix}{F}_{n}&{F}_{n + 1}&{F}_{n + 2}\end{bmatrix} = \begin{bmatrix}{F}_{1}&{F}_{2}&{F}_{3}\end{bmatrix} * {\begin{bmatrix}0&0&1\\1&0&0\\0&1&1\end{bmatrix}}^{n - 1}$$
看到这个矩阵递推公式就不用我多说什么了吧...
上AC代码:

//洛谷P1939
#include <stdio.h>
#include <string.h>
int tempa[3][3], tempb[3][3], M = 1000000007;
void matrixMultiply(int a[][3], int la, int ca, int b[][3], int lb, int cb)
{
    int i, j, k;
    for(i = 0; i < la; i++)
        for(j = 0; j < ca; j++)
            tempa[i][j] = a[i][j];
    for(i = 0; i < lb; i++)
        for(j = 0; j < cb; j++)
            tempb[i][j] = b[i][j];

    //清零
    for(i = 0; i < la; i++)
        for(j = 0; j < ca; j++)
            a[i][j] = 0;

    for(i = 0; i < la; i++)
        for(j = 0; j < cb; j++)
            for(k = 0; k < lb; k++){
                a[i][j] += (tempa[i][k] * (long long)(tempb[k][j])) % M;
                a[i][j] %= M;
            }
}
void fastPower(int base[][3], int p, int res[][3])
{
    while(p)
    {
        if(p & 1)   matrixMultiply(res, 1, 3, base, 3, 3);
        p >>= 1;
        matrixMultiply(base, 3, 3, base, 3, 3);
    }
}
int main(void)
{
    int t, n;
    scanf("%d", &t);
    while(t--)
    {
        int base[3][3] = {
            {0, 0, 1},
            {1, 0, 0},
            {0, 1, 1}
        };
        int res[3][3] = {1, 1, 1};
        scanf("%d", &n);
        if(n >= 4){
            fastPower(base, n - 1, res);
            printf("%d\n", res[0][0]);
        }else
            printf("1\n");
    }
    return 0;
}

本来还想再贴一题的,时间不允许。
大家一起加油!

声明:Trikker的小破站|版权所有,违者必究|如未注明,均为原创|本网站采用BY-NC-SA协议进行授权

转载:转载请注明原文链接 - 矩阵快速幂


分享一些学习成果