大一还是大二学线性代数的时候，感觉这玩意怎么这么无聊，这破矩阵能干啥？老师讲课跟催眠一样。后来学习OpenGL的时候才感受到矩阵的魅力，原来3D世界的平移、旋转、缩放都可以用一个小小的四阶矩阵的运算来完成，AMAZING 啊！！！(比老师上课有意思多了)
矩阵快速幂也是打ACM期间学到的，思想上和整数的快速幂一样，操作替换成了矩阵运算而已。
没错，课堂上的时间都用来睡觉、刷手机了
求解一些递推公式的第n项的时候，通过递推公式构造转移矩阵，并用矩阵快速幂可以快速得到第n项的值。特别对于n很大的时候不能循环迭代，可用矩阵快速幂解决。
比如斐波那契数列： $F_0=0,F_1=1,F_n=F_n-1+F_n-2,(n \ge2,n \in N^*)$
0 1 1 2 3 5 8…
如果要求第n项，n比较小的话，可以直接循环迭代出来。
如果n比较大，第一亿个是多少，用循环就太慢了，矩阵快速幂更快更适合。
根据矩阵乘法性质，构造转移矩阵
$$
\begin{bmatrix}
1&1\\
1&0\end{bmatrix} \times\begin{bmatrix}
F_{n-1}\\
F_{n-2}
\end{bmatrix} = \begin{bmatrix}
F_{n}\\
F_{n-1}
\end{bmatrix}
$$
那么
$$
\begin{bmatrix}
1&1\\
1&0\end{bmatrix}^{n-1} \times\begin{bmatrix}
F_{1}\\
F_{0}
\end{bmatrix} = \begin{bmatrix}
F_{n}\\
F_{n-1}
\end{bmatrix}
$$
这样问题就变成了求转移矩阵的幂。这个转移矩阵一定构造成方阵。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
|
#include <stdio.h>
#include <string.h>
#define MATRIX_SIZE 2
const int mod = 1e9+7;
struct Matrix //构造一个方阵
{
int data[MATRIX_SIZE][MATRIX_SIZE];
Matrix(){
memset(data, 0, sizeof(data));
for(int i=0; i<MATRIX_SIZE; i++){
data[i][i]=1; // 初始化为单位矩阵
}
}
Matrix operator + (Matrix o)const{ //(a+b)%mod
Matrix re;
for(int i=0; i<MATRIX_SIZE; i++){
for(int j=0; j<MATRIX_SIZE; j++){
re.data[i][j] = (this->data[i][j] + o.data[i][j]) % mod;
}
}
return re;
}
Matrix operator * (Matrix o)const{ //(a*b)%mod
Matrix re;
for(int i=0; i<MATRIX_SIZE; i++){
for(int j=0; j<MATRIX_SIZE; j++){
re.data[i][j] = 0;
for(int k=0; k<MATRIX_SIZE; k++)
re.data[i][j] = (re.data[i][j] + ((this->data[i][k] * o.data[k][j]) % mod)) % mod;
}
}
return re;
}
Matrix operator ^ (int n)const{ // (a^n)%mod
Matrix re,base;
base = *this;
while(n){
if(n&1)
re = re * base;
n>>=1;
base = base * base;
}
return re;
}
Matrix Psum(int n)const{ //(a+a^2+a^3.....+a^n)%mod
Matrix a,ans,pre;
int m;
a = *this;
if(n==1) return a;
m = n/2;
pre = a.Psum(m); // a^[1,n/2] 相加
ans = pre + (pre * (a^m)); // ans = [1,n/2] + a^(n/2)*[1,n/2]
if(n&1)
ans = ans + (a^n); //n为奇数时候a^n会漏掉，补上
return ans;
}
void out(){
for(int i=0; i<MATRIX_SIZE; i++){
for(int j=0; j<MATRIX_SIZE; j++){
printf("%d ",data[i][j]);
}printf("\n");
}
}
};
int main(){
Matrix ma,swa;
ma.data[0][0]=1;
ma.data[0][1]=1;
ma.data[1][0]=1;
ma.data[1][1]=0;
int n;
ma.out();
while(scanf("%d" , &n) , n>=1){
swa = ma^(n-1);
//swa.out();
printf("is %d\n", swa.data[0][0]*1 + swa.data[0][1]*0 );
}
return 0;
}
|
hdu6470
Farmer John有n头奶牛.
某天奶牛想要数一数有多少头奶牛,以一种特殊的方式:
第一头奶牛为1号,第二头奶牛为2号,第三头奶牛之后,假如当前奶牛是第n头,那么他的编号就是2倍的第n-2头奶牛的编号加上第n-1头奶牛的编号再加上自己当前的n的三次方为自己的编号.
现在Farmer John想知道,第n头奶牛的编号是多少,估计答案会很大,你只要输出答案对于123456789取模.
根据题意可以得到递推公式：
$$
F_n=\begin{cases}
1& n=1\\
2& n=2\\
2\times F_{n-2} + F_{n-1}+n^3& n\ge3
\end{cases}
$$
由于$n^3$的存在，这个递推公式并不是线性的，所以要想办法把$n^3$展开：
$$
n^3=(n-1+1)^3=C^0_3(n-1)^3+C^1_3(n-1)^2+C^2_3(n-1)^1+C^3_3(n-1)^0 \\=1\times(n-1)^3+3\times (n-1)^2+3\times(n-1)^1+1\times(n-1)^0
$$
那么转移矩阵$M$为：
$$
\begin{bmatrix}
0&1&0&0&0&0\\
2&1&C_3^0&C_3^1&C_3^2&C_3^3\\
0&0&C_3^0&C_3^1&C_3^2&C_3^3\\
0&0&0&C_2^0&C_2^1&C_2^2\\
0&0&0&0&C_1^0 &C_1^1\\
0&0&0&0&0&1
\end{bmatrix}\times
\begin{bmatrix}
F_{n-2}\\F_{n-1}\\(n-1)^3\\(n-1)^2\\(n-1)^1\\(n-1)^0
\end{bmatrix}
=
\begin{bmatrix}
F_{n-1}\\F_{n}\\(n)^3\\(n)^2\\(n)^1\\(n)^0
\end{bmatrix}
$$
取$n=3$：
$$
\begin{bmatrix}
F_{n-2}\\F_{n-1}\\(n-1)^3\\(n-1)^2\\(n-1)^1\\(n-1)^0
\end{bmatrix}=\begin{bmatrix}
F_{1}\\F_{2 }\\2^3\\ 2^2\\ 2^1\\ 2^0
\end{bmatrix}
$$
所以：
$$
F_n = M^{n-2 } \times \begin{bmatrix}
F_{1}\\F_{2 }\\2^3\\ 2^2\\ 2^1\\ 2^0
\end{bmatrix},n \ge 2
$$
AC代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
|
#include <stdio.h>
#include <string.h>
#include <iostream>
#define MATRIX_SIZE 6
using namespace std;
const int mod = 123456789;
struct Matrix //构造一个方阵
{
long long data[MATRIX_SIZE][MATRIX_SIZE];
Matrix(){
memset(data, 0, sizeof(data));
for(int i=0; i<MATRIX_SIZE; i++){
data[i][i]=1; // 初始化为单位矩阵
}
}
Matrix operator * (Matrix o)const{ //(a*b)%mod
Matrix re;
for(int i=0; i<MATRIX_SIZE; i++){
for(int j=0; j<MATRIX_SIZE; j++){
re.data[i][j] = 0;
for(int k=0; k<MATRIX_SIZE; k++)
re.data[i][j] = (re.data[i][j] + ((this->data[i][k] * o.data[k][j]) % mod)) % mod;
}
}
return re;
}
Matrix operator ^ (long long n)const{ // (a^n)%mod
Matrix re,base;
base = *this;
while(n){
if(n&1)
re = re * base;
n>>=1;
base = base * base;
}
return re;
}
};
Matrix ma,swa;
int main(){
int mm[6][6]={
{0,1,0,0,0,0},
{2,1,1,3,3,1},
{0,0,1,3,3,1},
{0,0,0,1,2,1},
{0,0,0,0,1,1},
{0,0,0,0,0,1}
};
for(int i=0;i<6;i++)
for(int j=0;j<6;j++)
ma.data[i][j]=mm[i][j];
long long t,a;
cin>>t;
while(t--){
cin>>a;
swa = ma^(a-2);
long long ac;
ac= swa.data[1][0]*1;
ac+=swa.data[1][1]*2;
ac+=swa.data[1][2]*8;
ac+=swa.data[1][3]*4;
ac+=swa.data[1][4]*2;
ac+=swa.data[1][5]*1;
cout<<ac%mod<<endl;
}
return 0;
}
|
这个方程怎么用矩阵计算，可以试着练练手。
$$
f(x)=
\begin{cases}
1& x=0 \\
1& x=1 \\
9\times f(x-1) + 9\times f(x-2) + 6& x>1
\end{cases}
$$