一些特殊的矩阵快速幂 hdu5950 hdu3369 hdu 3483

思想启发来自，

罗博士的根据递推公式构造系数矩阵用于快速幂

对于矩阵乘法和矩阵快速幂就不多重复了，网上很多博客都有讲解。主要来学习一下系数矩阵的构造

一开始，最一般的矩阵快速幂，要斐波那契数列Fn=F_n-1+F_n-2的第n项,想必都知道可以构造矩阵来转移

 其中，前面那个矩阵就叫做系数矩阵（我比较喜欢叫转移矩阵）

POJ3070 Fibonacci 可以试一试

#include<cstdio>
typedef long long ll;
const ll md=10000;
struct Mar{
    int r,c;
    ll a[10][10];
    Mar(){}
    Mar(int r,int c):r(r),c(c){
        for(int i=0;i<r;i++)
            for(int j=0;j<c;j++) a[i][j]=0;
    }
}A,T;
Mar mul(Mar A,Mar B){
    Mar ans(A.r,B.c);
    for(int i=0;i<A.r;i++)
        for(int j=0;j<B.c;j++)
            for(int k=0;k<A.c;k++){
                ans.a[i][j]+=A.a[i][k]*B.a[k][j]%md;
                if(ans.a[i][j]>=md) ans.a[i][j]-=md;
            }
    return ans;
}
Mar poww(Mar A,int b){
    Mar ans(A.r,A.c);
    for(int i=0;i<A.r;i++) ans.a[i][i]=1;
    while(b){
        if(b&1) ans=mul(ans,A);
        A=mul(A,A);
        b>>=1;
    }
    return ans;
}
void init(ll a,ll b){
    A=Mar(2,1);
    T=Mar(2,2);
    A.a[0][0]=b;A.a[1][0]=a;
    T.a[0][0]=T.a[0][1]=T.a[1][0]=1;
}
int main(){
    int t,n;
    ll a,b;
    while(~scanf("%d",&n)&&n!=-1){
        if(n==0) printf("0
");
        else{
            init(1,0);
            T=poww(T,n);
            A=mul(T,A);
            printf("%lld
",A.a[0][0]);
        }
    }
    return 0;
}

然后像类斐波那契数列Fn=a*Fn-1+b*Fn-2 ，也就变一下系数矩阵

而如果是要求斐波那契的前n项和，除去我们知道的第n+2项-1,那对于类斐波那契的前n项和，我们可以设S_n=S_n-1+F_{n，然后构造系数矩阵}

然后对于一些递推公式，转换一下前后几项的关系，也可以构造相应的系数矩阵来矩阵快速幂求解，比如求Σi^x的前n项和的话，设Fn=n^x,Sn=Sn-1+Fn

那么n^x可以写成(n-1+1)^x的形式，把n-1设为a，1设为b，然后就是二项式定理展开了，（二项式展开看右边→_→）

我们就可以得到一个n-1向n转移的递推过程，n^x=C^x_x(n-1)^x+C^x_x(n-1)^x-1+...+C⁰_x (n-1)⁰就可以构造系数矩阵

Recursive sequence HDU - 5950 

题意：F1=a,F2=b,Fn=Fn-1+2*Fn-2+n⁴,给出n求Fn

#include<cstdio>
typedef long long ll;
const ll md=2147493647;
struct Mar{
    int r,c;
    ll a[10][10];
    Mar(){}
    Mar(int r,int c):r(r),c(c){
        for(int i=0;i<r;i++)
            for(int j=0;j<c;j++) a[i][j]=0;
    }
}A,T;
Mar mul(Mar A,Mar B){
    Mar ans(A.r,B.c);
    for(int i=0;i<A.r;i++)
        for(int j=0;j<B.c;j++)
            for(int k=0;k<A.c;k++){
                ans.a[i][j]+=A.a[i][k]*B.a[k][j]%md;
                if(ans.a[i][j]>=md) ans.a[i][j]-=md;
            }
    return ans;
}
Mar poww(Mar A,int b){
    Mar ans(A.r,A.c);
    for(int i=0;i<A.r;i++) ans.a[i][i]=1;
    while(b){
        if(b&1) ans=mul(ans,A);
        A=mul(A,A);
        b>>=1;
    }
    return ans;
}
void init(ll a,ll b){
    A=Mar(7,1);
    T=Mar(7,7);
    A.a[0][0]=b;A.a[1][0]=a;
    for(int i=6,j=1;i>=2;i--,j<<=1) A.a[i][0]=j;
    T.a[0][0]=1;T.a[0][1]=2;T.a[1][0]=1;
    for(int i=6;i>=2;i--){
        T.a[i][6]=1;
        for(int j=5;j>=i;j--)
            T.a[i][j]=T.a[i+1][j]+T.a[i+1][j+1];
    }
    for(int i=2;i<7;i++) T.a[0][i]=T.a[2][i];
}
int main(){
    int t,n;
    ll a,b;
    scanf("%d",&t);
    while(t--){
        scanf("%d%lld%lld",&n,&a,&b);
        if(n==1) printf("%lld
",a%md);
        else if(n==2) printf("%lld
",b%md);
        else{
            init(a,b);
            T=poww(T,n-2);
            A=mul(T,A);
            printf("%lld
",A.a[0][0]);
        }
    }
    return 0;
}

那么(an+b)^x的前n项和呢，还是把n写成n-1+1的形式，然后(a(n-1)+1+b)^x二项式分解

(an+b)^x=C^x_x(n-1)^x*a^x*(a+b)⁰+C^x_x(n-1)^x-1*a^x-1*(a+b)¹+...+C⁰_x (n-1)⁰*a⁰*(a+b)^x,然后构造系数矩阵

RobotHDU - 3369 

题意：有个机器人第n天学n^k个单词，但星期六和星期天休息，不学单词，给出n，k,问机器一共学了多少个单词。

我们已经可以求n^k的前n项和了，那么接下来只需要求出星期六和星期天的再减去，即可。

那么，我们根据开始的星期几到下一个星期六距离的天数x,和下一个星期天距离的天数x+1，

那么要减去的部分就是，(7m1+x)^k 和(7m2+x+1)^k的前m项和m1=(n-x)/7,m2=(n-x-1)/7,且第一项分别是x^k和(x+1)^k

（转移矩阵懒得画了）

#include<cstdio>
typedef long long ll;
const ll md=1e9+7;
struct Mar{
    int r,c;
    ll a[25][25];
    Mar(){}
    Mar(int r,int c):r(r),c(c){
        for(int i=0;i<r;i++)
            for(int j=0;j<c;j++) a[i][j]=0;
    }
}A,T;
char s[25];
ll cf[25][25];
Mar mul(Mar A,Mar B){
    Mar ans(A.r,B.c);
    for(int i=0;i<A.r;i++)
        for(int j=0;j<B.c;j++)
            for(int k=0;k<A.c;k++){
                ans.a[i][j]+=A.a[i][k]*B.a[k][j]%md;
                if(ans.a[i][j]>=md) ans.a[i][j]-=md;
            }
    return ans;
}
Mar poww(Mar A,int b){
    Mar ans(A.r,A.c);
    for(int i=0;i<A.r;i++) ans.a[i][i]=1;
    while(b){
        if(b&1) ans=mul(ans,A);
        A=mul(A,A);
        b>>=1;
    }
    return ans;
}
void init(int n,ll f0,ll a,ll b){
    A=Mar(n,1);
    T=Mar(n,n);
    A.a[0][0]=f0,A.a[n-1][0]=1;
    for(int i=n-1;i>=1;i--){
        T.a[i][n-1]=1;
        for(int j=n-2;j>=i;j--){
            T.a[i][j]=T.a[i+1][j]+T.a[i+1][j+1];
            if(T.a[i][j]>=md) T.a[i][j]-=md;
        }
    }
    T.a[0][0]=1;
    for(int i=1;i<n;i++){
        T.a[0][i]=T.a[1][i];
        T.a[0][i]*=cf[a][n-1-i];
        if(T.a[0][i]>=md) T.a[0][i]%=md;
        T.a[0][i]*=cf[a+b][i-1];
        if(T.a[0][i]>=md) T.a[0][i]%=md;
    }
}
ll solve(int n,int m,int x){
    ll ans1,ans2,ans3;
    init(m+2,0,1,0);T=poww(T,n);
    A=mul(T,A);ans1=A.a[0][0];
    init(m+2,cf[x][m],7,x);T=poww(T,(n-x)/7);
    A=mul(T,A);ans2=A.a[0][0];
    init(m+2,cf[x+1][m],7,x+1);T=poww(T,(n-x-1)/7);
    A=mul(T,A);ans3=A.a[0][0];
    return ((((ans1-ans2)%md+md)%md-ans3)%md+md)%md;
}
int main(){
    for(int i=0;i<=20;i++){
        cf[i][0]=1;
        for(int j=1;j<=15;j++){
            cf[i][j]=cf[i][j-1]*i;
            if(cf[i][j]>=md) cf[i][j]%=md;
        }
    }
    int t=1,T,n,m,x;
    scanf("%d",&T);
    while(t<=T){
        scanf("%s",s);
        scanf("%d%d",&n,&m);
        if(s[0]=='M') x=6;
        else if(s[0]=='T'&&s[1]=='u') x=5;
        else if(s[0]=='W') x=4;
        else if(s[0]=='T') x=3;
        else if(s[0]=='F') x=2;
        else if(s[0]=='S'&&s[1]=='a') x=1;
        else x=0;
        ll ans=solve(n,m,x);
        printf("Case %d: %lld
",t++,ans);
    }
    return 0;
}

最后，如果是要求n^xxⁿ的前n项和呢，一样把n写成n-1+1，然后(n-1+1)^xxⁿ二项式分解

n^xxⁿ=C^x_x(n-1)^x*x^n-1+1+C^x_x(n-1)^x-1*x^n-1+1+...+C⁰_x (n-1)⁰*x^n-1+1,转移矩阵就有了

A Very Simple ProblemHDU - 3483 

求n^xxⁿ的前n项和的模板题

#include<cstdio>
typedef long long ll;
const int N=111;
struct Mar{
    int r,c;
    ll a[N][N];
    Mar(){}
    Mar(int r,int c):r(r),c(c){
        for(int i=0;i<r;i++)
            for(int j=0;j<c;j++) a[i][j]=0;
    }
}A,T;
ll md        ;
Mar mul(Mar A,Mar B){
    Mar ans(A.r,B.c);
    for(int i=0;i<A.r;i++)
        for(int j=0;j<B.c;j++)
            for(int k=0;k<A.c;k++){
                ans.a[i][j]+=A.a[i][k]*B.a[k][j]%md;
                if(ans.a[i][j]>=md) ans.a[i][j]-=md;
            }
    return ans;
}
Mar poww(Mar A,int b){
    Mar ans(A.r,A.c);
    for(int i=0;i<A.r;i++) ans.a[i][i]=1;
    while(b){
        if(b&1) ans=mul(ans,A);
        A=mul(A,A);
        b>>=1;
    }
    return ans;
}
void init(int n,int x){
    A=Mar(n,1);
    T=Mar(n,n);
    A.a[n-1][0]=1;
    T.a[0][0]=1;
    for(int i=n-1;i>=1;i--){
        T.a[i][n-1]=1;
        for(int j=n-2;j>=i;j--){
            T.a[i][j]=T.a[i+1][j]+T.a[i+1][j+1];
            if(T.a[i][j]>=md) T.a[i][j]-=md;
        }
    }
    for(int i=n-1;i>=1;i--)
        for(int j=n-1;j>=i;j--){
            T.a[i][j]*=x;
            if(T.a[i][j]>=md) T.a[i][j]%=md;
        }
    for(int i=1;i<n;i++) T.a[0][i]=T.a[1][i];
}
int main(){
    int n,x;
    while(~scanf("%d%d%lld",&n,&x,&md)){
        if(n<0&&x<0&&md<0) break;
        init(x+2,x);
        T=poww(T,n);
        A=mul(T,A);
        printf("%lld
",A.a[0][0]);
    }
    return 0;
}