快速幂算法及其在动态规划中的应用(矩阵幂)

魔术师的徒弟

已于 2022-01-29 21:51:05 修改

阅读量994

点赞数 1

分类专栏：算法专栏 LeetCode刷题文章标签：矩阵算法动态规划快速幂 C++

于 2022-01-29 18:29:37 首次发布

本文链接：https://blog.csdn.net/CS_COPy/article/details/122746335

版权

LeetCode刷题同时被 2 个专栏收录

64 篇文章 4 订阅

订阅专栏

算法专栏

27 篇文章 2 订阅

订阅专栏

文章目录

一、快速幂算法引入
二、快速幂算法
三、快速幂算法在矩阵幂中的运算
四、LeetCode1220的快速幂做法
五、参考资料

一、快速幂算法引入

假如我们有一个需要求2^100的后三位的问题，我们可以利用%运算的性质：
$b)\%p = (a\%p * b\%p)\%p$
这意味着我们如果想对 $2^{100}$ 直接取后三位，也就是取余1000，可以对它的每个因子取余1000,然后把这些因子乘起来再取余1000，得到的结果和直接对 $2^{100}$ 取余1000是一样的。

常规思路是循环100次，模拟乘100次2的过程。

但是这样思路写的代码经过时间测试用时有点长，我们能否优化这个O(N)的算法呢？

#include <iostream>
#include <ctime>
using namespace std;
long long normalPower(long long base, long long power)
{
    int res = 1;
    for (int i = 1; i <= power ; ++i)
    {
        res *= base;
        res %= 1000;
    }
    return res;
}

int main()
{
    clock_t start = clock();
    cout << normalPower(2, 1000000000) << endl;
    clock_t end = clock();
    cout << "The cost of calculate is " << double(end - start) / CLOCKS_PER_SEC << endl;
}

二、快速幂算法

核心思想就是如果幂次是偶数，那么我们可以让底数变成原来的平方，如果幂次是奇数，那么我们可以分解出一个底数，让它与返回结果相乘，然后幂次就是偶数了，通敌，让底数变为原来的平方，当幂次分到0的时候，就意味着我们不能再分解了，结束循环。

long long QuickPower(long long base, long long power)
{
    int res = 1;
    while (power > 0)
    {
        if (power % 2 == 0)
        {
            /*偶数情况 则底数平方 幂次变为原来的一半*/
            base = base * base % 1000;
            power /= 2;
        }
        else
        {
            /*指数是奇数情况 则拿出一个底数 使得指数变成偶数 然后就可以按照偶数那样操作*/
            res = res * base % 1000;
            --power;
            base = base * base % 1000;
            power /= 2;
        }
        /*直到幂次被削为0次幂时，循环结束，0次幂的时候就无法再进行削指数的操作了*/
    }
    return res;
}

优化一下重复部分，发现只有奇数有额外操作，可以简化为：

long long QuickPower(long long base, long long power)
{
    int res = 1;
    while (power > 0)
    {
        if (power % 2 == 1)
        {
            res = res * base % 1000;
            --power;
        }
        base = base * base % 1000;
        power /= 2;
        /*直到幂次被削为0次幂时，循环结束，0次幂的时候就无法再进行削指数的操作了*/
    }
    return res;
}

再优化一下，可以把power % 2 == 1优化为：power & 1 == 1，然后再优化成power & 1,然后power /= 2优化为power >>= 1

long long QuickPower(long long base, long long power)
{
    int res = 1;
    while (power > 0)
    {
        if (power & 1)
        {
            /*幂次为奇数时 需要把一个底数提取出来 然后让power - 1变为偶数*/
            res = res * base % 1000;
            --power;
        }
        /*偶数情况则直接让底数平方 指数除2*/
        base = base * base % 1000;
        power >>= 1;
        /*直到幂次被削为0次幂时，循环结束，0次幂的时候就无法再进行削指数的操作了*/
    }
    return res;
}

时间复杂度分析：快速幂算法的时间复杂度显然不是O(n)，而是O(logn)，算是很大的优化了。

三、快速幂算法在矩阵幂中的运算

以斐波那契数列为例，考虑行向量：(f(n),f(n-1)),由递推公式f(n) = f(n - 1) + f(n - 2),有：
$-1),f(n-2)]\begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}$
假设f(0)、f(1)已知，那么：
$[f(1),f(0)]\begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}^{n - 1}$
假设 $f (0) = 0, f (1) = 1$ ，我们的快速幂算法写法如下：

class Solution {
public:
    int Fibonacci(int n) 
    {
        if (n == 0)
        {
            return 0;
        }
        if (n == 1 || n == 2)
        {
            return 1;
        }
        int power = n - 1;
        vector<vector<long long>> base = {{1, 1}, {1, 0}};
        vector<vector<long long>> ret = {{1, 0}, {0, 1}};
        while (power > 0)
        {
            if (power & 1)
            {
                ret = mutiply(ret, base);
                --power;
            }
            base = mutiply(base, base);
            power >>= 1;
        }
        return ret[0][0];
        /*这里因为(f(1),f(0))是(1,0),那么f(n)就等于矩阵左上角元素*/
    }
    vector<vector<long long>> mutiply(const vector<vector<long long>>& a, 
                                      const vector<vector<long long>>& b)
    {
        vector<vector<long long>> c(2, vector<long long>(2));
        for (int i = 0; i < 2; ++i)
        {
            for (int j = 0; j < 2; ++j)
            {
                c[i][j] = a[i][0] * b[0][j] + a[i][1] * b[1][j];
            }
        }
        return c;
    }
};

四、LeetCode1220的快速幂做法

在那个题目中，我们得到了如下结果：

若定义向量 $\vec{f(i)}=(f(i,0),f(i,1),f(i,2),f(i,3),f(i,4))^{T}$ ，考虑矩阵知识，状态转移方程可以写为：
$\vec{f(i)}=\begin{bmatrix} 0 & 1&1&0&1 \\ 1 & 0&1&0&0\\ 0&0&1&0&0\\ 0&0&1&1&0\\ \end{bmatrix}\vec{f(i-1)}=A\vec{f(i-1)}$
且 $\vec{f(1)}=(1,1,1,1,1)^{T}$ ,要得到 $\vec{f(n)}$ 只要让f(1)乘 $A^{n-1}$ 即可.求矩阵的幂有一些快一点的数学算法，比如快速幂算法可以使用。

因此我们实现一下这个快速幂算法：

class Solution {
public:
    long long mod = 1e9 + 7;
    int countVowelPermutation(int n) 
    {
        if (n == 1)
        {
            return 5;
        }
        vector<vector<long long>> ret = {{1,0,0,0,0}, {0,1,0,0,0}, {0,0,1,0,0}, {0,0,0,1,0}, {0,0,0,0,1}};
        vector<vector<long long>> base = {{0,1,1,0,1}, {1,0,1,0,0}, {0,1,0,1,0}, {0,0,1,0,0}, {0,0,1,1,0}};
        int power = n - 1;
        while (power > 0)
        {
            if (power & 1)
            {
                ret = mutiply(base, ret);
                --power;
            }
            base = mutiply(base, base);
            power >>= 1;
        }
        long long ans = 0;
        /*此处对应A^(n-1)*f(1) 就等于所有元素之和*/
        for (int i = 0; i < 5; ++i)
        {
            for (int j = 0; j < 5; ++j)
            {
                ans = (ans + ret[i][j]) % mod;
            }
        }
        return ans % mod;
    }
    vector<vector<long long>> mutiply(const vector<vector<long long>>& a, const vector<vector<long long>>& b)
    {
        vector<vector<long long>> c(5, vector<long long>(5));
        for (int i = 0; i < 5; ++i)
        {
            for (int j = 0; j < 5; ++j)
            {
                c[i][j] = ((a[i][0] * b[0][j])%mod 
                + (a[i][1] * b[1][j])%mod 
                + (a[i][2] * b[2][j])%mod 
                + (a[i][3] * b[3][j])%mod
                + (a[i][4] * b[4][j])%mod)%mod;
            }
        }
        return c;
    }
};