目录

矩阵快速幂

hdu6470,斐波那契,矩阵快速幂,递推

解决哪些问题

求解一些递推公式的第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;
}