复习笔记之矩阵快速幂（不定时更新）

　　事情要从很久很久以前说起，我记得那是第一节算法课，老师说了这样一个题目：某人下楼梯可以一次走两阶，也可以一次走一阶，现在他在第 k 阶楼梯，问题是他要走到第0阶有多少种不同的走法？显然，这是fibonacci数列，即：

$$ f(n) = \begin{cases} 1& \text{n=1}\\ 1& \text{n=2}\\ f(n-1) + f(n-2)& \text(n \geq 3) \end{cases} $$

　　这个问题并不复杂，我们甚至能求出它的通项公式：

$$ \frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2})^{n}- (\frac{1-\sqrt{5}}{2})^{n}] $$

　　不过，我们知道斐波那契数列的数字总是整数，这个通项公式并没有显式的表明这是一个整数，而且，即使利用这个通项公式求解 $ f(n) $ 依然要计算至少 $ n $次，反倒不如直接递推！

　　事实上我们能在 log 级的时间复杂度内解决这个问题。

　　首先，我们构造矩阵：

$$ \begin{equation} \left| \begin{array}{cc} 1 & 1 \\ 1 & 0 \end{array} \right| \end{equation} $$

$$ \begin{equation} \left| \begin{array}{c} f(2) \\ f(1) \end{array} \right| \end{equation} $$

　　矩阵一和矩阵二相乘会得到一个新的矩阵：

$$ \begin{equation} \left| \begin{array}{c} f(2) + f(1) \\ f(2) \end{array} \right| \end{equation} $$

　　显然，矩阵三等价于：

$$ \begin{equation} \left| \begin{array}{c} f(3) \\ f(2) \end{array} \right| \end{equation} $$

　　如果把得到的矩阵三再和矩阵一相乘呢？我们会得到含有f(4)和f(3)的矩阵，依此类推，我们总能得到这样的矩阵：

$$ \begin{equation} \left| \begin{array}{c} f(n) \\ f(n-1) \end{array} \right| \end{equation} $$

　　记这个矩阵为$ M_{n} $ ，矩阵一为$ M_{0} $，矩阵二为$ M_{2} $则有：

$$ M_{n} = M_{0}*M_{n-1} = M_{0}*M_{0}*M_{n-2} = ... = M_{0}^{n-2}*M_{2} $$

　　作为一个"伪ACMer"我还是知道快速幂取模这个东西的：


　　 新得到的式子和快速幂取模能解决的问题竟如此相似，事实上是完全一样的，如果你了解快速幂取模的过程的话，唔---还是简要的说一下吧！


　　咦><整数和矩阵的区别很大么？把快速幂取模代码中的乘号重载成矩阵的乘法不就ok了吗？没错，就是这样！！！现在，已经成功的解决了这个问题！！

　　顺带说一下昨天BestCoder Round 80中的一个题目：

　　

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cstring>
#include <algorithm>
using namespace std ;
#define LL long long
#define rep(i,n) for (int i = 1 ; i <= n ; ++ i)
const int n = 3 ;
struct Matrix
{
LL mat[n<<1][n<<1] ;
}e;

void InitMatrix()
{
rep(i,n) rep(j,n) e.mat[i][j] = (i == j) ;
}

Matrix mul(Matrix a,Matrix b,LL Mod)
{
Matrix c ;
rep(i,n) {
rep(j,n) {
c.mat[i][j] = 0 ;
rep(k,n) {
c.mat[i][j] += (a.mat[i][k]*b.mat[k][j])%Mod ;
while (c.mat[i][j] >= Mod) c.mat[i][j] -= Mod ;
}
}
}
return c ;
}

Matrix Quick_Mod(Matrix a,LL k,LL Mod)
{
Matrix c = e ;
while (k) {
if (k&1) c = mul(c,a,Mod) ;
a = mul(a,a,Mod) ;
k = k/2 ;
}
return c ;
}

LL Quick_Mod(LL a,LL k,LL Mod)
{
LL ans = 1 ;
a %= Mod ;
while (k) {
if (k&1) ans = ans*a%Mod ;
a = a*a%Mod ;
k = k>>1 ;
}
return ans ;
}

int main()
{
//    freopen("in.txt","r",stdin) ;
LL N , A , B , C , Mod ;
int T ;
Matrix c ;
InitMatrix() ;
scanf("%d",&T) ;
while (T --) {
scanf("%I64d%I64d%I64d%I64d%I64d",&N,&A,&B,&C,&Mod) ;
if (A%Mod == 0) {
if (N == 1) puts("1") ;
else puts("0") ;
continue ;
}
if (N < 3) {
if (N == 1) printf("1\n") ;
else printf("%I64d\n",Quick_Mod(A,B,Mod)) ;
continue ;
}
memset(c.mat,0,sizeof(c.mat)) ;
c.mat[1][1] = C ;
c.mat[1][2] = c.mat[1][3] = c.mat[2][1] = c.mat[3][3] = 1 ;
c = Quick_Mod(c,N-2,Mod-1) ;
LL k = c.mat[1][1] + c.mat[1][3] ;
k %= (Mod-1) ;
k = k*B%(Mod - 1) ;
printf("%I64d\n",Quick_Mod(A,k,Mod)) ;
}
return 0 ;
}

HDU 5667

待学习链接：

　　matrix67的 十个利用矩阵乘法解决的经典题目

~end~