算法初步：快速乘，快速幂，矩阵快速幂

时间：2017-06-01 21:38:17 阅读：260 评论：0 收藏：0 [点我收藏+]

原创 by zoe.zhang

　　在刷题的时候遇到了问题，就是当循环或者递推的次数非常大的情况下获取一定结果，这个时候如果作普通运算，那么很容易就超时了，而且有时候结果也大得超范围了，即使是long long类型的也放不下，然后给了提示说是运用快速幂的思想。所以这里对快速幂做了一点思考和探讨。

1、快速乘，快速幂，矩阵快速幂三者的关系

技术分享

　　不管是快速乘，还是快速幂算法，实际上都包含了分解问题的思想在里面，将O（n）的复杂度降到O（lgn）。学习的时候，一般学习快速幂算法，再由此推广去解决矩阵快速幂问题。相对来说，矩阵快速幂应用更为广泛，在解决问题上更具有实用性，通常用来解决大数据的深度递推关系问题。（因为数据类型范围问题，所求的结果通常是取模之后的结果）矩阵快速幂与数值快速幂的关系就是将数值的相乘推广到矩阵之间的相乘，其本质思想是一致的。。

　　使用矩阵快速幂有两个步骤：

　　（1）从问题中提炼出递推转移方阵。（一般是方阵，方便计算和幂乘，这一步是解决问题的所在，因为一般来说矩阵快速幂可以套用程序模板）。

　　（2）采用快速幂的思想求解结果。

2、快速幂

　　那么什么是快速幂？字面意思，以较少的步骤、更有效率的方法获取到幂乘的结果。在这里可以通过两种角度来看快速幂乘法，两种角度其实是共通的，建议好好理解第一种角度。

　　角度一：幂指数的二进制化。

　　举个简单例子：a^9 = a ^ (1001)₂ = a^(2³*1+2²*0+2¹*0+2⁰*1) = (a^2³)*(a^2⁰) ;

　　快速幂的实际操作效果是减少了乘法的效果，且下一个乘数可以由上一个乘数的结果退出。实际代码如下所示，取模操作消耗比较大，采用位运算和&操作可以更好得提高程序效率。

long long quick_pow(int factor, int n)
{
    long long ans = 1;
    long long base = factor;
    while (n) 
    {
        if ((n & 1) != 0) ans *= base;  //奇偶判断
        base *= base;
        n >>= 1; //b/2  n>>=1 移位操作 注意 >>=
    }
    return ans;
}

　　这里用到了long long的类型，因为幂运算随随便便算起来得到的数都是很大的，有时候即使用了long long 型的数据还是放不下，所以通常会有一个取模的操作。而且通常快速幂的题目与取模是紧密相关的，当然关于取模的操作我还有一点知识点想捋一捋，所以放在下一节来说。

　　角度二：这个递推关系看起来比较明朗，相对第一种角度要好理解一点，本质上是层层分解，将很大的n值不断分解，并减少了乘法的次数。

技术分享

　　主要代码如下所示：

long long quick_pow2(int factor, int n)
{
    if (n == 0)
        return 1;
    if (n % 2 == 1)
        return quick_pow2(factor, n / 2)*quick_pow2(factor, n / 2) * factor;
    else
        return quick_pow2(factor, n / 2)*quick_pow2(factor, n / 2);
}

　　总的来说，递归的效率不如非递归的效率。所以建议使用第一种方法。

3、快速幂取模

　　好了，下面到了取模时间。实际上基本上90%的快速幂的题目需要取模。这里有一篇文章讲得通俗易懂，老少皆宜（https://wenku.baidu.com/view/d65f294702768e9951e73883.html）建议读一读。

　　取模公式：

（a+b）mod n = (（a mod n）+（b mod n）)mod n
（a-b） mod n = （（a mod n）-（b mod n）+n）mod n
ab mod n = （a mod n）（b mod n）mod n

　　这里注意一下减法那里 a mod n 可能会小于 b mod n 所以结果要加上n。

　　我们重点关注第三个公式，这是快速幂取模的理论基础。注意 ab mod n 可能在int里面，但是 a mod n 与 b mod n的乘积可能超过int 所以通常要用int来保存中间结果。

　　举一个笔试时遇到过的与快速幂有点关系的小题目，来加深对这个知识点的理解：

　　7^2014的十位数字是：(7^2014)%100 = (49^507)%100 = 49%100*((49*49%100)^253)%100 = 49%100 *(2401%100)^253 = 49%100 = 49;

　　快速幂取模的代码如下所示：

#define mod 1000000007
long long quick_pow(int factor, int n)
{
    long long ans = 1;
    long long base = factor; //这里也可以取一次模
    while (n)
    {
        if ((n & 1) != 0) ans = (ans * base) % mod;  //取模
        n >>= 1; //n/2  n>>=1 注意=
        base *= base;
        base %= mod;  //取模
    }
    return ans;
}

4、矩阵快速幂

　　OK，进入矩阵快速幂。矩阵快速幂实际上与数值快速幂没有什么区别，思想上是相同的，但是需要自己处理矩阵的乘法，通常方法是重载一下乘法运算符，或者是自定义矩阵乘法函数。需要注意一下的是，在作数值的快速幂时，基数ans=1，而矩阵的快速幂时，基数为单位矩阵，也就是ans = E。此外，程序中也需要注意一下矩阵乘法的自定义，有一些小的优化的点。

　　所以对于矩阵快速幂的题目，基本上就是套用一套矩阵快速幂方法，也就是通常所说的矩阵快速幂模板，因此解决问题的关键，就是找出正确的递推矩阵，这才是正确得到答案的关键所在。

　　大概就是所说，技巧是可以学习和磨练的，但是透过现象看到问题的本质才是根本和最重要的，也是在处理算法题目时成功的先决条件。

　　扯远了，从一些例子来学习一下矩阵快速幂。

　　例1：超级斐波那契数列:求斐波那契数列第一亿项的值（取模）。

技术分享

　　矩阵快速幂模板基本上是大同小异的，有一些细节要值得注意一下：例如：（1）减少取模运算，因为取模相对来说还是比较耗时的，（2）矩阵中0的个数比较多时候，可以对矩阵进行小幅度的优化。

#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;

const int N = 2, M = 2, P = 2;  //矩阵1 N*P； 矩阵2 :P *M
const int MOD= 100000007;
struct mat
{
    long long data[N][N];
};

mat A = { 1, 1, 1, 0 };
mat I = { 1, 0, 0, 1 };

//矩阵乘法定义
mat multi(mat a, mat b)
{
    mat ans;
    for (int i = 0; i < N; i++)
        for (int j = 0; j < M; j++)
        {
            ans.data[i][j] = 0;  //ans.data[][];
            for (int k = 0; k < P; k++)
                ans.data[i][j] += (a.data[i][k] * b.data[k][j] % MOD);//取模
            ans.data[i][j] %= MOD; //取模
        }
    return ans;
}

//矩阵快速幂
mat power(mat a, int n)  //a^K
{
    mat ans = I,base = a;
    while (n)
    {
        if (n & 1) ans = multi(ans, base);

        n >>= 1;//移位
        base = multi(base, base);
    }
    return ans;
}

int main()
{
    int n;
    cin >> n;
    mat ans = power(A, n - 1);  //n-1 注意一下 阶数
    cout << ans.data[0][0];
    system("pause");
    return 0;
}

　　下面是一个比较完善的快速幂模板，也是求解斐波那契数列的，但是更改矩阵可以套用在多个题目中，自己需要好好体会一下。

/快速幂模板
#include <iostream>
#include <cstdio>
using namespace std;

const int N = 2; //定义方阵的大小
const int MOD = 100000007;

struct mat{       //定义矩阵结构体；也可以定义一个类 加入成员变量等
    int data[N][N];
};

//注意快速幂里用作乘法的一般是方阵 所以简单起见
//矩阵A：N*N；矩阵B：N*N 所定义的结构体mat 也是一个方阵
//优化的写法
mat mul_1(mat A, mat B) //定义矩阵乘法
{
    mat ret;
    memset(ret.data, 0, sizeof(ret.data));  //初始化为0
    for (int i = 0; i < N;++i)
    {
        for (int k = 0; k < N; ++k)
        {
            if (A.data[i][k])  //为0
            {
                for (int j = 0; j < N; ++j)
                {
                    ret.data[i][j] += A.data[i][k] * B.data[k][j];
                    if (ret.data[i][j] >= MOD)   //减少一定的取模运算  取模操作耗时
                        ret.data[i][j] %= MOD;  
                }
            }
        }
    }
    return ret;
}

//未优化的写法
mat mul_2(mat A, mat B) //定义矩阵乘法
{
    mat ret;
    memset(ret.data, 0, sizeof(ret.data));  //初始化为0
    for (int i = 0; i < N; ++i)
    {
        for (int j = 0; j < N; ++j)
        {
                for (int k = 0; k < N; k++)
                {
                    ret.data[i][j] += A.data[i][k] * B.data[k][j];
                    //ret.data[i][j] += (A.data[i][k] * B.data[k][j])%MOD;
                    if (ret.data[i][j] >= MOD)   
                        ret.data[i][j] %= MOD;
                }
        }
    }
    return ret;
}

// 快速幂取模
mat power(mat A, int k)
{
    mat E, ans, base;
    
    if (k == 1) return A;

    memset(E.data, 0, sizeof(E.data));
    for (int i = 0; i < N; i++)
        E.data[i][i] = 1;

    if (k == 0) return E;
    
    ans = E; base = A;
    while (k)
    {
        if (k & 1) ans = mul_1(ans,base);
        base = mul_1(base, base);
        k >>= 1;
    }
    return ans;

}

int main()
{
    mat A = { 1, 1, 1, 0 };
    int k;
    while (cin>>k)
    {
        if (k == -1)break;
        mat ret = power(A, k);
        int ans = ret.data[0][1] % MOD;
        printf("%d\n", ans);
    }
    return 0;
}