Matrix Power Series

给出矩阵A，求矩阵\(A+A^2+...+A^k\)各个元素\(mod\ yyb\)的值，\(n\leq 30,k\leq 10^9,yyb\leq 10^4\)。

解

法一：分治

显然是数列题，故数列最浅显的减法是分治，寻找其中重复计算的部分，故可以

\[A+A^2+...+A^{k/2}+A^{k/2+1}+...+A^k=
\]

\[A+A^2+...+A^{k/2}+A^{k/2}(A+...+A^{k/2})+(k\&1)A^k=
\]

\[(A^{k/2}+1)(A+...+A^{k/2})+(k\&1)A^k
\]

对式子主体进行分治，其他的部分快速幂即可。

参考代码：

#include <iostream>

#include <cstdio>

#include <cstring>

#define il inline

#define ri register

using namespace std;

int n,yyb;

struct matrix{

    int jz[30][30];

    il void clear(){

        memset(jz,0,sizeof(jz));

    }

    il void unit(){

        clear();ri int i;

        for(i=0;i<n;++i)jz[i][i]|=true;

    }

    il void read(){

        ri int i,j;

        for(i=0;i<n;++i)

            for(j=0;j<n;++j)

                scanf("%d",&jz[i][j]);

    }

    il void print(){

        ri int i,j;

        for(i=0;i<n;++i,putchar('\n'))

            for(j=0;j<n;++j)

                printf("%d ",jz[i][j]);

        putchar('\n');

    }

    il matrix operator*(matrix x){

        matrix y;y.clear();ri int i,j,k;

        for(i=0;i<n;++i)

            for(j=0;j<n;y.jz[i][j]%=yyb,++j)

                for(k=0;k<n;++k)

                    y.jz[i][j]+=jz[i][k]*x.jz[k][j]%yyb;

        return y;

    }

    il matrix operator+(matrix x){

        matrix y;ri int i,j;

        for(i=0;i<n;++i)

            for(j=0;j<n;++j)

                y.jz[i][j]=(jz[i][j]+x.jz[i][j])%yyb;

        return y;

    }template<class free>

    il matrix operator^(free y){

        matrix x(*this),ans;ans.unit();

        while(y){

            if(y&1)ans=ans*x;

            x=x*x,y>>=1;

        }return ans;

    }

}state,unit,zero;

il matrix efen(int);

int main(){

    int k,i;

    scanf("%d%d%d",&n,&k,&yyb);

    unit.unit(),state.read();

    efen(k).print();

    return 0;

}

il matrix efen(int y){

    if(!y)return unit;

    if(y==1)return state;

    return ((state^(y>>1))+unit)*efen(y>>1)

        +((y&1)?(state^y):zero);

}

法二：矩阵快速幂

显然数列的题目，经常会存在递推方程，于是矩阵快速幂会在其中大有用武之地，于是设\(f_i=A+A^2+...+A^i\),不难有递推方程\(f_i=f_{i-1}+A^i\)，于是我们可以同时转移\(f_i,A^i\),故状态矩阵为

\[\begin{bmatrix}A_i&f_{i-1}\end{bmatrix}
\]

转移矩阵为

\[\begin{bmatrix}A&I\\0&I\end{bmatrix}
\]

以此矩阵中套矩阵仿照套路即可。

参考代码：

#include <iostream>

#include <cstdio>

#include <cstring>

#define il inline

#define ri register

using namespace std;

int n,yyb;

struct matrix1{

    int jz[30][30];

    il void read(){

        ri int i,j;

        for(i=0;i<n;++i)

            for(j=0;j<n;++j)

                scanf("%d",&jz[i][j]);

    }

    il void print(){

        ri int i,j;

        for(i=0;i<n;++i,putchar('\n'))

            for(j=0;j<n;++j)

                printf("%d ",jz[i][j]);

        putchar('\n');

    }

    il void clear(){

        memset(jz,0,sizeof(jz));

    }

    il void unit(){

        clear();ri int i;

        for(i=0;i<n;++i)jz[i][i]|=true;

    }

    il matrix1 operator*(matrix1 x){

        matrix1 y;y.clear();

        ri int i,j,k;

        for(i=0;i<n;++i)

            for(j=0;j<n;y.jz[i][j]%=yyb,++j)

                for(k=0;k<n;++k)

                    y.jz[i][j]+=jz[i][k]*x.jz[k][j]%yyb;

        return y;

    }

    il matrix1 operator+(matrix1 x){

        matrix1 y;ri int i,j;

        for(i=0;i<n;++i)

            for(j=0;j<n;++j)

                y.jz[i][j]=(jz[i][j]+x.jz[i][j])%yyb;

        return y;

    }template<class free>

    il matrix1 operator^(free y){

        matrix1 ans,x(*this);ans.unit();

        while(y){

            if(y&1)ans=ans*x;

            x=x*x,y>>=1;

        }return ans;

    }

}A;

struct matrix2{

    matrix1 jz[2][2];

    il void clear(){

        memset(jz,0,sizeof(jz));

    }

    il void unit(){

        clear();ri int i;

        for(i=0;i<2;++i)jz[i][i].unit();

    }

    il matrix2 operator*(matrix2 x){

        matrix2 y;y.clear();

        ri int i,j,k;

        for(i=0;i<2;++i)

            for(j=0;j<2;++j)

                for(k=0;k<2;++k)

                    y.jz[i][j]=y.jz[i][j]+jz[i][k]*x.jz[k][j];

        return y;

    }template<class free>

    il matrix2 operator^(free y){

        matrix2 ans,x(*this);ans.unit();

        while(y){

            if(y&1)ans=ans*x;

            x=x*x,y>>=1;

        }return ans;

    }

}state,tran;

int main(){

    int k;

    scanf("%d%d%d",&n,&k,&yyb);

    A.read(),state.jz[0][0]=A;

    tran.jz[0][0]=A,tran.jz[0][1].unit();

    tran.jz[1][1].unit(),state=state*(tran^k);

    state.jz[0][1].print();

    return 0;

}

法三：倍增前缀和优化

数列求和可以利用倍增优化，主要思想是维护\(A^{2^i}\)，这个显然可通过\(A^{2^{i+1}}=(A^{2^i})^2\)来暴力递推，再维护一个倍增的数列，\(f_i=A^1+A^2+...+A^{2^i}\)，不难得知其转移方程为\(f_i=f_{i-1}A^{2^{i-1}}+f_{i-1}\)，而这个的维护也可以暴力维护，于是对于我们的求和，以\(A+A^2+...+A^{15}\)为例，有

\[A+A^2+...+A^{15}=f_3+A^9+...+A^{15}=
\]

\[f_3+A^8(A+...+A^7)=f_3+A^8(f_2+A^5+A^6+A^7)=
\]

\[f_3+A^8[f_2+A^4(A+A^2+A^3)]=
\]

\[f_3+A^8[f_2+A^4(f_1+A^2A^1)]=f_3+A^8[f_2+A^4(f_1+A^2f_0)]
\]

于是按照这个先维护好一段倍增前缀和，再对后式提公因式，继续利用倍增前缀和优化，不停地继续，即可得到答案，当然你也可以递归处理，以下参考程序用的是非递归的方式。

参考代码：

#include <iostream>

#include <cstdio>

#include <cstring>

#define il inline

#define ri register

using namespace std;

int n,yyb;

struct matrix{

    int jz[30][30];

    il void clear(){

        memset(jz,0,sizeof(jz));

    }

    il void unit(){

        clear();ri int i;

        for(i=0;i<n;++i)jz[i][i]|=true;

    }

    il void read(){

        ri int i,j;

        for(i=0;i<n;++i)

            for(j=0;j<n;++j)

                scanf("%d",&jz[i][j]);

    }

    il void print(){

        ri int i,j;

        for(i=0;i<n;putchar('\n'),++i)

            for(j=0;j<n;++j)

                printf("%d ",jz[i][j]);

        putchar('\n');

    }

    il matrix operator*(matrix x){

        matrix y;ri int i,j,k;y.clear();

        for(i=0;i<n;++i)

            for(j=0;j<n;y.jz[i][j]%=yyb,++j)

                for(k=0;k<n;++k)

                    y.jz[i][j]+=jz[i][k]*x.jz[k][j]%yyb;

        return y;

    }

    il matrix operator+(matrix x){

        matrix y;ri int i,j;

        for(i=0;i<n;++i)

            for(j=0;j<n;++j)

                y.jz[i][j]=(jz[i][j]+x.jz[i][j])%yyb;

        return y;

    }

}A,p[31],sum[31],ans,jilu;

int main(){

    int k,i;

    scanf("%d%d%d",&n,&k,&yyb);

    A.read(),p[0]=A,jilu.unit(),sum[0]=A;

    for(i=1;i<=30;++i)p[i]=p[i-1]*p[i-1];

    for(i=1;i<=30;++i)sum[i]=sum[i-1]*p[i-1]+sum[i-1];

    for(i=30;i>=0;--i)

        if(k>>i&1){

            ans=ans+sum[i]*jilu;

            jilu=jilu*p[i];

        }ans.print();

    return 0;

}

总上所诉，不难得知数列题目的一般求法，矩阵快速幂和分治，而特殊地对于数列前缀和的问题，我们可以利用倍增前缀和优化，但无论如何，数列问题绝不只一条道路通往罗马。

巴特西

Matrix Power Series

解

最新文章

热门文章