摘要 初识矩阵快速幂已是很久之前的事情了,我的忘性有那么好,现在自然是忘得一干二净啦,简要写一下矩阵快速幂相关,以备后用.

  事情要从很久很久以前说起,我记得那是第一节算法课,老师说了这样一个题目:某人下楼梯可以一次走两阶,也可以一次走一阶,现在他在第 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"我还是知道快速幂取模这个东西的:

快速幂取模能在log级的时间复杂度内解决 $ A^{n} mod p $的问题, A , p , n 均为正整数

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

  代码是这样的:

 

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

 

  思路是这样的:如果 n 是奇数,那么$ answer = A^{n-1}*A mod p  $,如果 n 是偶数,那么$ answer = (A^{2})^{(n/2)} mod p $。这样,就会把时间复杂度降到log级别

 

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

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

  

 

  对,就是这个,读完题目我就知道要构造矩阵了,瞎写了一个多小时,都没能正确的构造矩阵,也因此才写了这篇随笔。

  分析 我们给它的递推式换一种写法:

$$ f_{n} = \begin{cases} f_{2}^{0}-->0& \text{n=1}\\ f_{2}^{1}-->1& \text{n=2}\\ f_{2}*f_{n-1}^{c}*f_{n-2}-->(1+c*g_{n-1}+g_{n-2})& \text(n \geq 3) \end{cases} $$

  这样写的意思就是都表示为$ f_{2} $的某次方,最后得到$ g_{n} = 1 + c*g_{n-1} + g_{n-2} $且有了$ g_{1} = 0, g_{2} = 1 $

  于是构造矩阵:

$$ \left| \begin{array}{ccc} c & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array} \right| $$

$$ \left| \begin{array}{c} g_{2} \\ g_{1} \\ 1 \end{array} \right| $$

  可以发现,相乘一次就能得到$ g_{3} $。同样,到这里这个问题就解决了,下面是代码。

我的代码在HDoj上提交可以 AC ,但是直觉告诉我应该有点什么问题==希望路过的同学能帮忙指出!

 

 1 #include <iostream>
 2 #include <stdio.h>
 3 #include <stdlib.h>
 4 #include <cstring>
 5 #include <algorithm>
 6 using namespace std ;
 7 #define LL long long
 8 #define rep(i,n) for (int i = 1 ; i <= n ; ++ i)
 9 const int n = 3 ;
10 struct Matrix
11 {
12     LL mat[n<<1][n<<1] ;
13 }e;
14 
15 void InitMatrix()
16 {
17     rep(i,n) rep(j,n) e.mat[i][j] = (i == j) ;
18 }
19 
20 Matrix mul(Matrix a,Matrix b,LL Mod)
21 {
22     Matrix c ;
23     rep(i,n) {
24         rep(j,n) {
25             c.mat[i][j] = 0 ;
26             rep(k,n) {
27                 c.mat[i][j] += (a.mat[i][k]*b.mat[k][j])%Mod ;
28                 while (c.mat[i][j] >= Mod) c.mat[i][j] -= Mod ;
29             }
30         }
31     }
32     return c ;
33 }
34 
35 Matrix Quick_Mod(Matrix a,LL k,LL Mod)
36 {
37     Matrix c = e ;
38     while (k) {
39         if (k&1) c = mul(c,a,Mod) ;
40         a = mul(a,a,Mod) ;
41         k = k/2 ;
42     }
43     return c ;
44 }
45 
46 LL Quick_Mod(LL a,LL k,LL Mod)
47 {
48     LL ans = 1 ;
49     a %= Mod ;
50     while (k) {
51         if (k&1) ans = ans*a%Mod ;
52         a = a*a%Mod ;
53         k = k>>1 ;
54     }
55     return ans ;
56 }
57 
58 
59 int main()
60 {
61 //    freopen("in.txt","r",stdin) ;
62     LL N , A , B , C , Mod ;
63     int T ;
64     Matrix c ;
65     InitMatrix() ;
66     scanf("%d",&T) ;
67     while (T --) {
68         scanf("%I64d%I64d%I64d%I64d%I64d",&N,&A,&B,&C,&Mod) ;
69         if (A%Mod == 0) {
70             if (N == 1) puts("1") ;
71             else puts("0") ;
72             continue ;
73         }
74         if (N < 3) {
75             if (N == 1) printf("1\n") ;
76             else printf("%I64d\n",Quick_Mod(A,B,Mod)) ;
77             continue ;
78         }
79         memset(c.mat,0,sizeof(c.mat)) ;
80         c.mat[1][1] = C ;
81         c.mat[1][2] = c.mat[1][3] = c.mat[2][1] = c.mat[3][3] = 1 ;
82         c = Quick_Mod(c,N-2,Mod-1) ;
83         LL k = c.mat[1][1] + c.mat[1][3] ;
84         k %= (Mod-1) ;
85         k = k*B%(Mod - 1) ;
86         printf("%I64d\n",Quick_Mod(A,k,Mod)) ;
87     }
88     return 0 ;
89 }
HDU 5667

 

 待学习链接:

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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

~end~