【数论】矩阵快速幂&&CODE[VS] 1732

矩阵，快速幂以及矩阵快速幂终于到了最后一个重头戏

---

矩阵快速幂

（矩阵与快速幂的无限基情）
矩阵快速幂基于矩乘和快速幂的一种算法
它的一个作用是加速递推，例如快速求出斐波那契数列

关于矩阵乘法：http://blog.youkuaiyun.com/loi_peacefuldog/article/details/52783928

关于快速幂：http://blog.youkuaiyun.com/loi_peacefuldog/article/details/52783942

而矩阵快速幂，其实是一个很水的东西
就是将快速幂的数值替换为一个矩阵，原理完全相同，但还是有些小地方需要注意

首先是一种特殊的矩阵——单位矩阵，这是一种除了主对角线上数值为1，矩阵其余各点，数值均为0的正方（行列数相同）矩阵

１　０　０　０　０　０
０　１　０　０　０　０
０　０　１　０　０　０
０　０　０　１　０　０
０　０　０　０　１　０
０　０　０　０　０　１　

例如这一个矩阵，单位矩阵的大小是任意的，根据具体题意来定，例如在求斐波那契数列第n项时，行列数当等于为2，不然不满足矩阵乘法的先决条件

那么这个单位矩阵有什么用处呢，我们可以把它当做实数里面的1，我们知道任意一个实数乘以1等于它本身，任意一个实数乘以它的倒数（1/该实数）等于1，矩阵也具有相似的性质，即任意一个矩阵乘以单位矩阵都等于原矩阵本身

所以我们那可以知道：
单位矩阵* 单位矩阵 * ….* 单位矩阵 == 单位矩阵
单位矩阵 => 单位矩阵 * 单位矩阵 => 单位矩阵的^2 * 单位矩阵^2 => 单位矩阵^4 * 单位矩阵^4…….
单位矩阵^n = 单位矩阵^(2^t1)* 单位矩阵^(2^t2) * …… * 单位矩阵^(2^tn)

友情链接（快速幂）：http://blog.youkuaiyun.com/loi_peacefuldog/article/details/52783942

这样一来便可以用单位矩阵来替代快速幂中进行快速连乘的xx,上面的博文解释了这样做的原理，这里就不再赘述（其实就是懒得打字了）

另外一个需要注意就是递推公式的求解和何如将递推公式转化为矩阵
（递推公式这需要自己根据题意推导，感觉没法讲…..)

拿斐波那契数列为例
我们都知道
斐波那契数列的递推公式为
Ｆ［ｎ］　＝　｛　
　　　　　　　　　０　（ｎ　＝　０）
　　　　　　　　　１　（ｎ　＝　１）
　　　　　　　　　Ｆ［n－１］＋Ｆ［n－２］　（ｎ　＞　１）
　　　　　　　｝；

那么怎么将这个公式转化为矩阵呢

这里用到斐波那契数列的一个性质：即前n项fibonacci的平方和等于第fib[n]*fib[n+1]。

我们先来想F[ ｎ ] = Ｆ［ｎ－１］＋Ｆ［ｎ－２］　（ｎ　＞　１）时矩阵应当如何构造

咳咳，显然我们很容易想到一个行数为1，列数为2的矩阵

| F［ｎ－１］　　Ｆ［ｎ－２］|

然后显然：
已知F[ｎ] = Ｆ［ｎ－１］＋Ｆ［ｎ－２］
则可以写成F[ｎ]　= 　Ｆ［ｎ－１］＊１＋Ｆ［ｎ－２］＊１
　　　　　　　　＝　（Ｆ［ｎ－１］＋Ｆ［ｎ－２］）＊１

因为F［１］ = F［２］ = 1,所以 又可以写成F[ｎ] = Ｆ［ｎ－１］＊Ｆ［１］＋Ｆ［ｎ－２］＊Ｆ［２］

？？？你是不是已经发现了什么？？？
对呀这个对应关系是符合矩阵乘法运算法则的

所以我们便又可以写出：

｜Ｆ［ｎ］　　＝　｜１　１｜＊｜Ｆ［ｎ－１］　
｜Ｆ［ｎ－１］　　｜１　０｜　｜Ｆ［ｎ－２］

用数学表达式即是：
①F［ｎ］　＝　Ｆ［ｎ－１］＋Ｆ［ｎ－２］
②F［ｎ－１］　＝ Ｆ［ｎ－１］＋０

但是我们想算某个Ｆ［ｎ］但是不知道Ｆ［ｎ－１］与Ｆ［ｎ－２］的值怎么办呢？

因为本身F[ ｎ ] = Ｆ［ｎ－１］＋Ｆ［ｎ－２］ 便是一个递推公式
所以我们可以将将这个式子化到最初状态即：
　　　　　　　　　　　　　　　　　　　　Ｆ［３］＝Ｆ［１］＋Ｆ［２］
　　　　　　　　　　　　　　　　　　　　　（２　＝　１　＋　１）
化为矩阵则为：

｜Ｆ［３］　＝　｜１　１｜＊｜Ｆ［２］　
｜Ｆ［２］　　　｜１　０｜　｜Ｆ［１］

将Ｆ［３］和Ｆ［２］换为F[ ｎ ] 和　Ｆ［ｎ－１］
由已知的递推关系可得:

｜Ｆ［ｎ］　　＝　｜１　１｜＾（ｎ－２）＊｜Ｆ［２］　
｜Ｆ［ｎ－１］　　｜１　０｜　　　　　　　｜Ｆ［１］

接下来直接递推求得Ｆ［ｎ］即可了
再利用快速幂加速递推，使得计算步数大大减少，是不是非常exciting

在代码实现上最好用一个结构体将数组都封装起来
这样做不仅能使代码更加好看，也更加易读更容易调试和挑出bug
并且封装之后代码量大大降低，写起来更简单
（一开始也是直接打的二维数组代码很乱直到看见一篇神犇的博客）
（http://www.cnblogs.com/vongang/archive/2012/04/01/2429015.html orz）

附上模板水题代码：CODE[VS] 1732(POJ 3070)

#include <cstdio>
#include <cstring> 
#include <iostream>
#include <algorithm>
#define MOD 1000000007

using namespace std;
typedef long long LL; 

LL n;

struct osu{
    LL m[2][2];
}ans, csj;

osu jc(osu a, osu b){
    osu tmp;
    for(LL i = 0; i <= 1;i++){
        for(LL j = 0; j <= 1;j++){
            tmp.m[i][j] = 0;
            for(LL k = 0; k <= 1;k++)
                tmp.m[i][j] = (tmp.m[i][j]% MOD + (a.m[i][k]%MOD * b.m[k][j]%MOD)%MOD) % MOD;
//(a * b) % p = (a % p * b % p) % p,(a + b) % p = (a % p + b % p) % p
        }
    }
    return tmp;
}
LL pow(LL n){
    csj.m[0][0] = 1;
    csj.m[0][1] = 1;
    csj.m[1][0] = 1;
    csj.m[1][1] = 0;
    for (LL i = 0;i <= 1;i++)  
        for (LL j = 0;j <= 1;j++){
            ans.m[i][j] = 0; 
        } 
    for (LL i = 0;i <= 1;i++)  
        ans.m[i][i] = 1;
    while(n){
        if(n & 1){
            ans = jc(ans, csj);
        }
        n >>= 1;
        csj = jc(csj, csj);
    }
    return ans.m[0][1];
}

int main(){
    while(scanf("%lld",&n) != EOF){   
        printf("%lld\n",pow(n));
    }

return 0;
}

---

THE END

蛤蛤，终于完了 准备开始搞几个题，然后继续搞数论
（预定再搞2天差不多了，明天还要考试求RP++）

By Peacefuldoge
http://blog.youkuaiyun.com/loi_peacefuldog