[模板]矩阵的快速幂

 

 原理：

$\left( {\begin{array}{*{20}{c}}
{{F_{n + 2}}}\\
{{F_{n + 1}}}
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
1&1\\
1&0
\end{array}} \right)\left( \begin{array}{l}
{F_{{\rm{n + }}1}}\\
{F_n}
\end{array} \right)$

 记这个矩阵为A，则有

$\left( {\begin{array}{*{20}{c}}
{{F_{n + 1}}}\\
{{F_n}}
\end{array}} \right) = {A^n}\left( \begin{array}{l}
{F_1}\\
{F_0}
\end{array} \right)$

 

注意这种初始化方式，matrix为只有一个二维数组的结构体

struct matrix
{
    long long m[4][4];
};
matrix ans = {  1,0,0,0,
                0,1,0,0,
                0,0,1,0,
                0,0,0,1,};
matrix base = { 1,1,1,0,
                1,0,0,0,
                0,1,0,0,
                1,0,0,1,};

 

 

第一种要注意，界限的确定，要保证符合所有矩阵。

 

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;
struct mat{
    ll m[2][2];
};
ll n=10;
ll mod=10000;
mat mul(mat &A,mat &B){
    mat C={0};
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            for(int k=0;k<2;k++){
                C.m[i][j]=(C.m[i][j]+A.m[i][k]*B.m[k][j])%mod;
            }
        }
    }
    return C;
}
mat pow(mat A,ll n){
    mat B={0};
    B.m[0][0]=B.m[1][1]=1;
    while(n>0){
        if(n&1) B=mul(B,A);
        A=mul(A,A);
        n>>=1;
    }
    return B;
}
int main(){
    mat A={0};
    A.m[0][0]=A.m[1][0]=A.m[0][1]=1;
    A=pow(A, n);
    printf("%lld\n",A.m[1][0]);
}

 

 

 

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef vector<ll> vec;
typedef vector<vec> mat;
ll n=10;
const int M=10000;
mat mul(mat &A,mat &B){
    mat C(A.size(),vec(B[0].size()));
    for(int i=0;i<(int)A.size();i++){
        for(int j=0;j<(int)B[0].size();j++){
            for(int k=0;k<(int)B.size();k++){
                C[i][j]=(C[i][j]+A[i][k]*B[k][j])%M;
            }
        }
    }
    return C;
} 

mat pow(mat A,ll n){
    mat B(A.size(),vec(A.size()));
    for(int i=0;i<(int)A.size();i++){
        B[i][i]=1;//单位矩阵 
    }
    while(n>0){
        if(n&1) B=mul(B,A);
        A=mul(A,A);
        n>>=1;
    }
    return B;
}
int main(){
    mat A(2,vec(2));
    A[0][0]=1;A[0][1]=1;
    A[1][0]=1;A[1][1]=0;
    A=pow(A,n);
    printf("%lld",A[1][0]);//为什么输出这个，按照递推公式得出 
     
}

 

 

http://acm.xidian.edu.cn/problem.php?id=1145

 

再介绍一题http://acm.hit.edu.cn/hoj/problem/view?id=2255

是哈工大的在线oj上的一个题目，一个类似于斐波那契的题目，题目会给出a,b,p,q,s,e, 其中f(0)=a,f(1)=b,当n>=2时 f(n)=P*f(n-1)+q*f((n-2) 求它组成的数列从第s项起一直加到第e项的和是多少，这题就不能推f(n)的矩阵乘法的递推式了，得推出他的前n项和s(n)的递推式,当然如果会推斐波那契的矩阵形式的递推形式的想必这个就不难了

推理过程如下：

$F\left( N \right) = S\left( N \right) - S\left( {N - 1} \right);$

又有$F\left( N \right) = P*F\left( {N - 1} \right) + Q*F\left( {N - 2} \right)$

所以 $S\left( N \right) - S\left( {N - 1} \right) = P*\left( {S\left( {N - 1} \right) - S\left( {N - 2} \right)} \right) + Q*\left( {S\left( {N - 2} \right) - S\left( {N - 3} \right)} \right)$

$S\left( N \right) = \left( {P + 1} \right)*S\left( {N - 1} \right) + \left( {Q - P} \right)*S\left( {N - 2} \right) - Q*S\left( {N - 3} \right)$

推到这里你应该知道怎么把它化成矩阵乘法的形式了吧？注意到右边有三个S(x)项所以矩阵递推式左边也要有三项写出来就是

化简一下

所以第s项到第e项的和就是S(E)-S(s-1) 注意是S(s-1)

 

解题方式与上面类似

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;
ll mod=1000000009;
struct mat{
    ll m[4][4];
};
ll l,r,res1,res2;
mat A={ 2,0,0,-1,
        1,0,0,0,
        0,1,0,0,
        0,0,1,0};
mat mul(mat &A,mat &B){
    mat C={0};
    for(int i=0;i<4;i++){
        for(int j=0;j<4;j++){
            for(int k=0;k<4;k++){
                C.m[i][j]=(C.m[i][j]+A.m[i][k]*B.m[k][j]%mod+mod)%mod;
            }
        }
    }
    return C;
}

mat mod_pow(mat A,ll n){
    mat B={0};
    B.m[0][0]=B.m[1][1]=B.m[2][2]=B.m[3][3]=1;
    while(n>0){
        if(n&1) B=mul(B,A);
        A=mul(A,A);
        n>>=1;
    }
    return B;
}

ll res(mat &a){
    return (a.m[3][0]*6+a.m[3][1]*3+a.m[3][2]*2+a.m[3][3])%mod;
}

int main(){
    while(scanf("%lld%lld",&l,&r)!=EOF){
        if(l==0) res1=0;
        else{
            mat l2=mod_pow(A,l-1);
            res1=res(l2);
        }
        if(r==0) res2=1;
        else{
        mat rr=mod_pow(A,r);
            res2=res(rr);
        }
        ll ans=(res2-res1+mod)%mod;
        printf("%lld\n",ans);
    }
}