【浮*光】#矩阵乘法# 矩阵优化的学习与练习

T：【p4838】p哥破解密码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<string>
#include<cmath>
using namespace std;
typedef long long ll;

/*【p4838】p哥破解密码 //煞笔题我调了一个小时 /哭哭...
定义一个串合法，当且仅当串只由A和B构成，且没有连续的3个A。
密码就是长度为N的合法字符串数量对19260817取模的结果。*/

/*【分析】f[i][0]=f[i-1][0]+f[i-1][1]+f[i-1][2];
f[i][1]=f[i-1][0],f[i][2]=f[i-1][1]; f[1][0]=f[1][1]=1; */

/* 
得到dp矩阵：1 1 0
           1 0 1
           1 0 0 
那么答案就是：矩阵(1,1,0)*dp矩阵^(n-1)得到的a.m[1][1]+a.m[1][2]+a.m[1][3]。*/

//【矩阵快速幂】给定n*n的矩阵A，求A^k。

void reads(ll &x){ //读入优化（正负整数）
    ll fx=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fx=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=x*10+s-'0';s=getchar();}
    x*=fx; //正负号
}

const ll mod=19260817;

struct Mat{ ll m[9][9]; }a,dp; //a是原始矩阵

Mat mat_mul(Mat x,Mat y){ //矩阵乘：x*y=c
    Mat c; memset(c.m,0,sizeof(c.m)); //记得要随时清空
    for(ll i=1;i<=3;i++) for(ll j=1;j<=3;j++) for(ll k=1;k<=3;k++)
        c.m[i][j]=(c.m[i][j]%mod+x.m[i][k]*y.m[k][j]%mod)%mod; return c; }

Mat mat_pow(Mat x,ll y) //矩阵快速幂：ans*x^y
 {  Mat ans; memset(ans.m,0,sizeof(ans.m)); //记得要随时清空
    for(ll p=1;p<=3;p++) ans.m[p][p]=1; //初始化为单位矩阵（一定不能忘记233）
    while(y){ if(y&1) ans=mat_mul(ans,x); x=mat_mul(x,x); y>>=1; } return ans; }

void init()
 { memset(a.m,0,sizeof(a.m)); a.m[1][1]=1,a.m[1][2]=1,a.m[1][3]=0;
   memset(dp.m,0,sizeof(dp.m)); dp.m[1][1]=1,dp.m[1][2]=1,dp.m[2][1]=1,dp.m[2][3]=1,dp.m[3][1]=1; }

int main(){ 
    ll T,n; reads(T); while(T--){
        reads(n); init(); Mat tmp=mat_pow(dp,n-1); a=mat_mul(a,tmp);
        cout<<(a.m[1][1]+a.m[1][2]+a.m[1][3])%mod<<endl; }
}

T：【p4910】帕秋莉的手环

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<string>
#include<cmath>
using namespace std;
typedef long long ll;

/*【p4910】帕秋莉的手环 */

/*【分析】f[i][0]=f[i-1][0]+f[i-1][1](=f[i-2][0]); 
f[i][1]=f[i-1][0]; f[1][0]=f[1][1]=1; 发现这其实是一个斐波那契数列。
所以答案为：求斐波那契数列 第n项 + 第n-1项的两倍 。*/

/* 得到dp矩阵：1 1 
              0 1 。*/


//  3 4  乘  递推矩阵  =  4 7 （手算也行）

//  解上述方程得递推矩阵为： 

//             0  1
//  3 4  乘            =  4 7
//             1  1

void reads(ll &x){ //读入优化（正负整数）
    ll fx=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fx=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=x*10+s-'0';s=getchar();}
    x*=fx; //正负号
}

const ll mod=1000000007;

struct Mat{ ll m[9][9]; }a,dp; //a是原始矩阵

Mat mat_mul(Mat x,Mat y){ //矩阵乘：x*y=c
    Mat c; memset(c.m,0,sizeof(c.m)); //记得要随时清空
    for(ll i=1;i<=2;i++) for(ll j=1;j<=2;j++) for(ll k=1;k<=2;k++)
        c.m[i][j]=(c.m[i][j]%mod+x.m[i][k]*y.m[k][j]%mod)%mod; return c; }

Mat mat_pow(Mat x,ll y) //矩阵快速幂：ans*x^y
 {  Mat ans; memset(ans.m,0,sizeof(ans.m)); //记得要随时清空
    for(ll p=1;p<=2;p++) ans.m[p][p]=1; //初始化为单位矩阵（一定不能忘记233）
    while(y){ if(y&1) ans=mat_mul(ans,x); x=mat_mul(x,x); y>>=1; } return ans; }

void init()
 { memset(dp.m,0,sizeof(dp.m)); dp.m[1][1]=0,dp.m[1][2]=dp.m[2][1]=dp.m[2][2]=1;
   memset(a.m,0,sizeof(a.m)); a.m[1][1]=2,a.m[1][2]=1,a.m[2][1]=a.m[2][2]=0; }

int main(){ 
    ll T,n; reads(T); while(T--){
        reads(n); init(); Mat tmp=mat_pow(dp,n-1); a=mat_mul(a,tmp);
        cout<<(a.m[1][2])%mod<<endl; } //这个位置就是 第n项 + 第n-1项的两倍
}

T：【sp2742】SumSums

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<string>
#include<cmath>
using namespace std;
typedef long long ll;

/*【sp2742】SumSums
长度为N的整数序列C，进行T次操作，每次把序列每一位修改为'其他各位'的和。 */

/*  每一对 原数 和 加密后的数 构成一个矩阵，即矩阵a：[ci,sum-ci]。
    则加密一次后的矩阵A为[sum-ci,(n-1)sum-(sum-ci)]。
    可以发现：所有数加密一次后总和变为原来的n-1倍。
    推出 a 乘 dp矩阵[[0,n-1],[1,n-2]] 可以得到 A。
    按照这个规律，加密T次和T+1次构成的矩阵为a*b^t。
    对于每个数，处理出n个矩阵a，就可以用快速幂解决a*bt，得到加密t次的数。*/

void reads(ll &x){ //读入优化（正负整数）
    ll fx=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fx=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=x*10+s-'0';s=getchar();}
    x*=fx; //正负号
}

const ll mod=98765431; ll c[50019],sum=0;

struct Mat{ ll m[9][9]; }a,dp;

Mat mat_mul(Mat x,Mat y){ //矩阵乘：x*y=c
    Mat c; memset(c.m,0,sizeof(c.m)); //记得要随时清空
    for(ll i=1;i<=2;i++) for(ll j=1;j<=2;j++) for(ll k=1;k<=2;k++)
        c.m[i][j]=(c.m[i][j]%mod+x.m[i][k]*y.m[k][j]%mod)%mod; return c; }

Mat mat_pow(Mat x,ll y) //矩阵快速幂：ans*x^y
 {  Mat ans; memset(ans.m,0,sizeof(ans.m)); //记得要随时清空
    for(ll p=1;p<=2;p++) ans.m[p][p]=1; //初始化为单位矩阵（一定不能忘记233）
    while(y){ if(y&1) ans=mat_mul(ans,x); x=mat_mul(x,x); y>>=1; } return ans; }

int main(){ 
    ll n,t; reads(n),reads(t); //T次操作 
    for(ll i=1;i<=n;i++) reads(c[i]),sum=(sum+c[i])%mod;
    dp.m[1][1]=0,dp.m[1][2]=n-1,dp.m[2][1]=1,dp.m[2][2]=n-2;
    dp=mat_pow(dp,t); //先求出dp矩阵的T次方
    for(ll i=1;i<=n;i++){ a.m[2][1]=a.m[2][2]=0;
        a.m[1][1]=c[i],a.m[1][2]=(sum-c[i]+mod)%mod;
        printf("%lld
",mat_mul(a,dp).m[1][1]);
    }
}

T：【p3702】序列计数

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<string>
#include<queue>
#include<map>
#include<vector>
#include<cmath>
using namespace std;
typedef long long ll;

//【p3702】序列计数
//符合要求：序列中的数都是不超过m的正整数，且这n个数的和是p的倍数。
//Alice还希望，这n个数中，至少有一个数是质数。有多少个序列满足她的要求。

//至少有1个质数的方案数=总方案数-不含质数的方案数
//可以先在O(m)的时间内把1~m的质数筛出来。
//然后考虑：从mod p=0到mod p=0，可以由mod p=a和mod p=p-a两个阶段组成。
//设f[i][j]表示从mod p=i到mod p=j的方案数，这是一个矩阵，自乘m次就能得到答案。

//对于每个数可以更新所有的f[i]（0≤i<n），不过这样会TLE。
//因为f都是循环出现的，只需要求f[0]即可，再据此推其余的f。

void reads(int &x){ //读入优化（正负整数）
    int fa=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fa=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+s-'0';s=getchar();}
    x*=fa; //正负号
}

const int mod=20170408;

bool vis_[20000019]; int primes[5000019],cnt;

struct Mat{ int n,m; ll num[110][110];
    Mat(){ n=m=0; memset(num,0,sizeof(num)); }
    Mat operator*(const Mat a)const{
      Mat ans; ans.n=n,ans.m=a.m; int i,j,k;
      for(i=1;i<=ans.n;i++) for(j=1;j<=ans.m;j++) for(k=1;k<=m;k++)
        ans.num[i][j]=(ans.num[i][j]+num[i][k]*a.num[k][j])%mod; return ans; } }A,B;

Mat mat_pow(Mat x,int y)
 { Mat ans; ans.n=x.n,ans.m=x.m; for(int i=1;i<=ans.n;i++) ans.num[i][i]=1;
   while(y){ if(y&1) ans=ans*x; x=x*x,y>>=1; } return ans; }

int main(){
    int n,m,p; reads(n),reads(m),reads(p); vis_[1]=1;
    for(int i=2;i<=m;i++){ //欧拉线性筛质数
        if(!vis_[i]) primes[++cnt]=i;
        for(int j=1;j<=cnt&&i*primes[j]<=m;j++)
         { vis_[i*primes[j]]=1; if(i%primes[j]==0) break; }
    } A.n=A.m=B.n=B.m=p; //设置为p*p的矩阵
    for(int i=1;i<=m;i++){ A.num[p][(i-1)%p+1]++; //先推p
        if(vis_[i]) B.num[p][(i-1)%p+1]++; } //所有数 vs 合数
    for(int i=p-1;i>=1;i--) for(int j=1;j<=p;j++) //线性逆推
        A.num[i][j]=A.num[i+1][j%p+1],B.num[i][j]=B.num[i+1][j%p+1];
    printf("%lld
",(mat_pow(A,n).num[p][p]-mat_pow(B,n).num[p][p]+mod)%mod);
}

T：【p2461】递归数列

// luogu-judger-enable-o2
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<string>
#include<cmath>
using namespace std;
typedef long long ll;

/*【p2461】递归数列
对于i<=k：ai=bi；对于i>k:ai=c1*ai-1+c2*ai-2+...+ck*ai-k。
给定自然数m<=n, 计算am+a(m+1)+a(m+2)+...+an, 并输出它%p的余数的值。*/

/* 对于矩阵[a1 a2 ... ak]，按照题目中的公式推出[a2 a3 ... ak+1]。
   将区间求和转化为前缀相减处理(注意此时不要改变dp矩阵a，需要特判t<=k)。*/
   
//真·调了一个下午，矩阵题就是这么...emm...

//最后是因为在 int main() 里也写了一个 ll p; ...一直re

void reads(ll &x){ //读入优化（正负整数）
    ll fx_=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fx_=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+s-'0';s=getchar();}
    x*=fx_; //正负号
}

ll n,p,sum[20];

struct Mat{ ll v[20][20]; //矩阵的数值v[][]
    Mat(){ memset(v,0,sizeof(v)); }
    Mat operator*(const Mat a)const{ Mat ans; ll i,j,k;
      for(i=1;i<=n;i++) for(j=1;j<=n;j++) for(k=1;k<=n;k++)
        ans.v[i][j]=(ans.v[i][j]+v[i][k]*a.v[k][j]%p)%p; return ans; } }A,B;

Mat mat_pow(Mat x,ll y)
 { Mat ans; for(ll i=1;i<=n;i++) ans.v[i][i]=1;
   while(y){ if(y&1) ans=ans*x; x=x*x,y>>=1; } return ans; }

ll cal(ll x){ if(x<n) return sum[x]; return (B*mat_pow(A,(x-n+1))).v[1][n]; }

int main(){
    reads(n); for(int i=1;i<=n;i++) reads(B.v[1][i]), //n+1行存前缀和
        B.v[1][n+1]=B.v[1][n+1]+B.v[1][i],sum[i]=sum[i-1]+B.v[1][i];
    for(int i=n;i>=1;i--) reads(A.v[i][n]),A.v[i][n+1]=A.v[i][n];
    for(int i=1;i<n;i++) A.v[i+1][i]=1; n++; A.v[n][n]=1;
    ll l,r; reads(l),reads(r),reads(p);
    for(int i=1;i<n;i++) sum[i]%=p;
    for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) 
        A.v[i][j]%=p,B.v[i][j]%=p;
    printf("%lld
",(cal(r)-cal(l-1)+p)%p);
}

T：【p2151】HH去散步

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<string>
#include<cmath>
using namespace std;
typedef long long ll;

/*【p2151】HH去散步
N个点M条边的无向图，A->B。不能立刻沿着上次边的反方向走。求距离为t方案数。*/

/*【分析】把每条边看作行列序号。则(i,j)表示：从一条边走到另一条边的方案数。
对于两条边i,j，如果不是刚走过来的边并且能够接上(i.y=j.x)，则Aij=1。*/

void reads(int &x){ //读入优化（正负整数）
    int fx_=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fx_=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+s-'0';s=getchar();}
    x*=fx_; //正负号
}

const int p=45989,N=159; int n,m,k,s,t,px[N],py[N];

struct Mat{ int v[N][N]; //矩阵的数值v[][]
    Mat(){ memset(v,0,sizeof(v)); }
    Mat operator*(const Mat a)const{ Mat ans; int i,j,k;
      for(i=0;i<=m;i++) for(j=0;j<=m;j++) for(k=0;k<=m;k++)
        ans.v[i][j]=(ans.v[i][j]+v[i][k]*a.v[k][j]%p)%p; return ans; } }a;

Mat mat_pow(Mat x,int y)
 { Mat ans; for(int i=0;i<=m;i++) ans.v[i][i]=1;
   while(y){ if(y&1) ans=ans*x; x=x*x,y>>=1; } return ans; }

int main(){
    reads(n),reads(m),reads(k),reads(s),reads(t);
    for(int i=1;i<=m;i++) reads(px[i<<1]),reads(py[i<<1]), 
        px[i<<1|1]=py[i<<1],py[i<<1|1]=px[i<<1]; //反向边
    m=m*2+1; for(int i=2;i<=m;i++){
        if(px[i]==s) a.v[0][i]=1; if(py[i]==t) a.v[i][1]=1;
        for(int j=2;j<=m;j++) //(i.y=j.x)，可以接上这个边
            if(py[i]==px[j]&&(i^j)!=1) a.v[i][j]=1; //不是同一条边
    } a=mat_pow(a,k+1); printf("%d
",a.v[0][1]);
}

T：【p2151】迷路

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<queue>
#include<string>
#include<cmath>
using namespace std;
typedef long long ll;

/*【p2151】迷路
n个点的有向图，每条边的权值都在[1,9]之间。求‘从1到n经过路径边权和恰好为t’的方案数模2009。*/

/*【分析】考虑拆点，将1个点拆成9个连成链，每次将‘入点i对应边权的点’与‘出点j’相连。*/

//https://www.cnblogs.com/GXZlegend/p/8302711.html

void reads(int &x){ //读入优化（正负整数）
    int fx_=1;x=0;char s=getchar();
    while(s<'0'||s>'9'){if(s=='-')fx_=-1;s=getchar();}
    while(s>='0'&&s<='9'){x=(x<<3)+(x<<1)+s-'0';s=getchar();}
    x*=fx_; //正负号
}

const int p=2009,N=109; int n,m,k;

struct Mat{ int v[N][N]; //矩阵的数值v[][]
    Mat(){ memset(v,0,sizeof(v)); }
    Mat operator*(const Mat a)const{ Mat ans; int i,j,k;
      for(i=0;i<=m;i++) for(j=0;j<=m;j++) for(k=0;k<=m;k++)
        ans.v[i][j]=(ans.v[i][j]+v[i][k]*a.v[k][j]%p)%p; return ans; } }A;

Mat mat_pow(Mat x,int y)
 { Mat ans; for(int i=0;i<=m;i++) ans.v[i][i]=1;
   while(y){ if(y&1) ans=ans*x; x=x*x,y>>=1; } return ans; }

int main(){
    reads(n),reads(k); m=n*9; char ss[19];
    for(int i=0;i<n;i++){ scanf("%s",ss);
        for(int j=0;j<8;j++) A.v[i*9+j][i*9+j+1]=1;
        for(int j=0;j<n;j++) //↑↑拆成的9个点连成一条链
            if(ss[j]!='0') A.v[i*9+ss[j]-'1'][j*9]=1; 
        //有向图(i,j)连边，选取入点i的第ss[j]-'1'个记录边权，出点在j*9位置
    } printf("%d
",mat_pow(A,k).v[0][n*9-9]); // 0 -> n-1
}

T：【p3990】超级跳马

#include <cstdio>
#include <cstring>
#include <algorithm>
#define mod 30011
using namespace std;

//【p3990】[Shoi2013]超级跳马 //矩阵乘法

int n; //见：https://www.cnblogs.com/GXZlegend/p/7816112.html

struct data{
    int v[105][105];
    data() {memset(v , 0 , sizeof(v));}
    int *operator[](int a) {return v[a];}
    data operator*(data &a){
        data ans; int i , j , k;
        for(i=1;i<=n;i++) for(j=1;j<=n;j++) for(k=1;k<=n;k++)
            ans[i][j]=(ans[i][j]+v[i][k]*a[k][j]) % mod; return ans;
    }
}I , A , B;

data pow(data x , int y){
    data ans;
    for(int i = 1 ; i <= n ; i ++ ) ans[i][i] = 1;
    while(y){
        if(y & 1) ans = ans * x;
        x = x * x , y >>= 1;
    } return ans;
}

int main(){
    int m , i; scanf("%d%d" , &n , &m);
    for(i = 1 ; i <= n ; i ++ ) I[i][i] = I[i + n][i] = I[i][i + n] = 1;
    for(i = 1 ; i < n ; i ++ ) I[i + 1][i] = I[i][i + 1] = 1;
    n <<= 1 , A = pow(I , m - 2) , B = A * I;
    printf("%d
" , (B[1][n >> 1] - A[1][n] + mod) % mod);
}

 

                                                ——时间划过风的轨迹，那个少年，还在等你