洛谷题单 Part 6.7.1 矩阵

应队友要求，开始学线性代数，具体路线是矩阵$\to$高斯消元$\to$线性基。为多项式做个准备

P3390 【模板】矩阵快速幂

题面
板子，用结构体写的，感觉有点丑，一会儿看看题解有没有写得好看的

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 110;
const ll mod=1e9+7;
struct node{ll a[N][N];int len;}sqr;
void sqr0(node &x){memset(x.a,0,sizeof x.a);x.len=sqr.len;
}
void sqr1(node &x){memset(x.a,0,sizeof x.a);x.len=sqr.len;for(int i=1;i<=x.len;i++)x.a[i][i]=1;
}
node operator*(node x, node b){node c;sqr0(c);for(int i=1;i<=x.len;i++){for(int j=1;j<=x.len;j++){for(int k=1;k<=x.len;k++)(c.a[i][j]+=x.a[i][k]*b.a[k][j]%mod)%=mod;}}return c;
}void qpow(node &x, ll y){node re;sqr1(re);while(y){if(y&1)re=re*x;x=x*x;y>>=1;}x=re;
}
ll k;
int main(){scanf("%d%lld",&sqr.len,&k);for(int i=1;i<=sqr.len;i++){for(int j=1;j<=sqr.len;j++)scanf("%lld",&sqr.a[i][j]);}qpow(sqr,k);for(int i=1;i<=sqr.len;i++){for(int j=1;j<=sqr.len;j++)printf("%lld ",sqr.a[i][j]);puts("");}
}

P1939 【模板】矩阵加速（数列）

题面
搞个方阵
$$A_3=\left[\begin{array}{ccc}{a}_3& a_2& a_1\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],X=\left[\begin{array}{ccc}1& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right],$$
则$$A_{3}X=\left[\begin{array}{ccc}{a}_4& a_3& a_2\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]=A_4,$$
因此对$X$进行矩阵快速幂即可。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5;
const ll mod=1e9+7;
struct node{ll a[N][N];}sqr,A;
void sqr0(node &x){memset(x.a,0,sizeof x.a);
}
void sqr1(node &x){memset(x.a,0,sizeof x.a);for(int i=1;i<=3;i++)x.a[i][i]=1;
}
node operator*(node x, node b){node c;sqr0(c);for(int i=1;i<=3;i++){for(int j=1;j<=3;j++){for(int k=1;k<=3;k++)(c.a[i][j]+=x.a[i][k]*b.a[k][j]%mod)%=mod;}}return c;
}void qpow(node &x, ll y){node re;sqr1(re);while(y){if(y&1)re=re*x;x=x*x;y>>=1;}x=re;
}
ll n,T;
int main(){cin>>T;while(T--){cin>>n;if(n<=3){puts("1");continue;}sqr0(sqr);sqr.a[1][1]=sqr.a[1][2]=sqr.a[2][3]=sqr.a[3][1]=1;sqr0(A);A.a[1][1]=A.a[1][2]=A.a[1][3]=1;qpow(sqr,n-3);A=A*sqr;cout<<A.a[1][1]<<endl;}}

P4783 【模板】矩阵求逆

题面
把一个矩阵通过行变换变为单位矩阵所需要的行变换操作，操作给一个单位矩阵，就可以得到其逆矩阵。故应用高斯消元即可。

#include<bits/stdc++.h>
#define N 1000
using namespace std;
const int mod=1e9+7;
inline void read(int &x){int s=0,w=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}while(ch>='0'&&ch<='9'){s=(s<<3)+(s<<1)+(ch&15);ch=getchar();}x=s*w;
}
int n,a[N][N];
int qpow(int x, int y){int re=1;while(y){if(y&1)re=1LL*re*x%mod;x=1LL*x*x%mod,y>>=1;}return re;
}
int main(){read(n);for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)read(a[i][j]);a[i][n+i]=1;}for(int i=1;i<=n;i++){int now=i;for(int j=i;j<=n;j++)if(a[now][i]<a[j][i])now=j;if(a[now][i]==0){puts("No Solution");return 0;}if(now!=i)swap(a[now],a[i]);for(int j=i+1;j<=n<<1;j++)a[i][j]=1LL*a[i][j]*qpow(a[i][i],mod-2)%mod;a[i][i]=1;for(int j=1;j<=n;j++){if(j==i)continue;int div=1LL*a[j][i]*qpow(a[i][i],mod-2)%mod;for(int k=i;k<=n<<1;k++)a[j][k]=(a[j][k]-1LL*a[i][k]*div%mod+mod)%mod;}}for(int i=1;i<=n;i++,puts(""))for(int j=1;j<=n;j++)printf("%d ",a[i][n+j]);}

P1962 斐波那契数列

题面
构造矩阵
$$A_2=\left[\begin{array}{cc}{f}_2& f_1\\ 0& 0\end{array}\right],X=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right],$$
则$$A_{2}X=\left[\begin{array}{cc}{f}_3& f_2\\ 0& 0\end{array}\right]=A_3,$$

#include<cstdio>
typedef long long ll;
const ll mod=ll(1e9+7);
struct node
{ll sqr[5][5];
}a;
node operator*(node a, node b)
{node c;c.sqr[1][1]=(a.sqr[1][1]*b.sqr[1][1]%mod+a.sqr[1][2]*b.sqr[2][1]%mod)%mod;c.sqr[1][2]=(a.sqr[1][1]*b.sqr[1][2]%mod+a.sqr[1][2]*b.sqr[2][2]%mod)%mod;c.sqr[2][1]=(a.sqr[2][1]*b.sqr[1][1]%mod+a.sqr[2][2]*b.sqr[2][1]%mod)%mod;c.sqr[2][2]=(a.sqr[2][1]*b.sqr[1][2]%mod+a.sqr[2][2]*b.sqr[2][2]%mod)%mod;return c;
}
ll n;
void quickpow(node &x, ll y)
{node rec;rec.sqr[1][1]=rec.sqr[2][2]=1,rec.sqr[1][2]=rec.sqr[2][1]=0;while(y){if(y&1)rec=rec*x;x=x*x,y>>=1;}x=rec;
}
int main()
{scanf("%lld",&n);if(n==0)return puts("0");a.sqr[1][1]=a.sqr[1][2]=a.sqr[2][1]=1,a.sqr[2][2]=0;quickpow(a,n-1);printf("%lld\n",a.sqr[1][1]);
}

P1349 广义斐波那契数列

题面
构造矩阵
$$A_2=\left[\begin{array}{cc}{f}_2& f_1\\ 0& 0\end{array}\right],X=\left[\begin{array}{cc}P& 1\\ Q& 0\end{array}\right],$$
则$$A_{2}X=\left[\begin{array}{cc}{f}_3& f_2\\ 0& 0\end{array}\right]=A_3,$$

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll mod;
struct node
{ll sqr[5][5];node(){memset(sqr,0,sizeof sqr);}
}a,b;
node operator*(node a, node b)
{node c;c.sqr[1][1]=(a.sqr[1][1]*b.sqr[1][1]%mod+a.sqr[1][2]*b.sqr[2][1]%mod)%mod;c.sqr[1][2]=(a.sqr[1][1]*b.sqr[1][2]%mod+a.sqr[1][2]*b.sqr[2][2]%mod)%mod;c.sqr[2][1]=(a.sqr[2][1]*b.sqr[1][1]%mod+a.sqr[2][2]*b.sqr[2][1]%mod)%mod;c.sqr[2][2]=(a.sqr[2][1]*b.sqr[1][2]%mod+a.sqr[2][2]*b.sqr[2][2]%mod)%mod;return c;
}
ll n;
void quickpow(node &x, ll y)
{node rec;rec.sqr[1][1]=rec.sqr[2][2]=1,rec.sqr[1][2]=rec.sqr[2][1]=0;while(y){if(y&1)rec=rec*x;x=x*x,y>>=1;}x=rec;
}
int main()
{scanf("%lld%lld%lld%lld%lld%lld",&a.sqr[1][1],&a.sqr[2][1],&b.sqr[1][2],&b.sqr[1][1],&n,&mod);if(n<=2)return printf("%lld\n",b.sqr[1][3-n]);a.sqr[1][2]=1,a.sqr[2][2]=0;quickpow(a,n-2);b=b*a;printf("%lld\n",b.sqr[1][1]);
}

P4000 斐波那契数列

题面
不是，这什么题都往题单里放啊，这是我$18$年外出集训堵了个论文费了两三天时间才切了的人生中第一道黑题，现在变成紫题了。
有一个性质是$f_n\mathrm{mod}p$有循环节，且循环节长度不会超过$6p$，这还有个名叫皮萨诺定理。所以我们考虑求出循环节的长度，然后用矩阵乘法求出结果。
引理：对于$f_n\mathrm{mod}p$的循环节$g(p)$有如下性质：
$1.p=p_i^{{\alpha }_i}$，即$p$为质数的幂时，$g(p)=g(p_i)\times p_i^{{\alpha }_i-1}$
$2.p=\prod p_i^{{\alpha }_i}$，即$p$为合数时，$g(p)=lcm(g(p_i^{{\alpha }_i})$
对于$g(p)$这么算，如果$5$是模$p$的二次剩余，那么循环节为$p-1$的因子，否则为$2p+2$的因子。
因为$p$不是特别大，直接取$p-1$和$2p+2$即可。
对于$p\le 5$就暴力算即可，$g(2)=3,g(3)=5,g(5)=20$

// luogu-judger-enable-o2
#include<cstdio>
#include<cmath>
#include<cstring>
using namespace std;
#define rg register
typedef long long ll;
char str[30000000];
ll n,p,mod,len,fac[100000],power[100000],faccnt,s;
struct node
{ll sqr[5][5];
}b;
node operator *(node a, node b)
{node xx;xx.sqr[1][1]=(a.sqr[1][1]*b.sqr[1][1]%p+a.sqr[1][2]*b.sqr[2][1]%p)%p;xx.sqr[1][2]=(a.sqr[1][1]*b.sqr[1][2]%p+a.sqr[1][2]*b.sqr[2][2]%p)%p;xx.sqr[2][2]=(a.sqr[2][1]*b.sqr[1][2]%p+a.sqr[2][2]*b.sqr[2][2]%p)%p;xx.sqr[2][1]=(a.sqr[2][1]*b.sqr[1][1]%p+a.sqr[2][2]*b.sqr[2][1]%p)%p;return xx;
}
node quickpow(node x, ll y)
{node rec;rec.sqr[1][1]=rec.sqr[2][2]=1;rec.sqr[1][2]=rec.sqr[2][1]=0;while(y){if(y%2==1)rec=rec*x;x=x*x;y/=2;}return rec;
}
ll gcd(ll a, ll b)
{if(b==0)return a;else return gcd(b,a%b);
}
ll lcm(ll a, ll b)
{return a*b/gcd(a,b);
}
ll get(ll k)
{ll now=k;for(rg ll i=2;i*i<=now;i++){if(now%i==0){faccnt++;fac[faccnt]=i;power[faccnt]=1;while(now%i==0){now/=i;power[faccnt]*=i;}}}for(rg ll i=1;i<=faccnt;i++)power[i]/=fac[i];if(now!=1){fac[++faccnt]=now;power[faccnt]=1;}for(rg ll i=1;i<=faccnt;i++){if(fac[i]==2)power[i]*=3;else if(fac[i]==3)power[i]*=5;else if(fac[i]==5)power[i]*=20;else if(fac[i]%5==1||fac[i]%5==4)power[i]*=fac[i]-1;else power[i]*=(fac[i]+1)<<1;}ll ans=power[1];for(rg ll i=1;i<=faccnt;i++)ans=lcm(ans,power[i]);return ans;
}
int main()
{scanf("%s%lld",str,&p);if(p==1){printf("0\n");return 0;}mod=get(p);len=strlen(str);for(rg ll i=0;i<len;i++)n=((n<<3)+(n<<1)+(str[i]&15))%mod;if(n==0){printf("0\n");return 0;}if(n==1||n==2){printf("1\n");return 0;}b.sqr[2][2]=0;b.sqr[1][1]=b.sqr[1][2]=b.sqr[2][1]=1;b=quickpow(b,n-1);printf("%lld\n",b.sqr[1][1]);
}

P3758 [TJOI2017] 可乐

题面

#include<bits/stdc++.h>
#define N 50
using namespace std;
const int mod=2017;
int t,n,m;
struct node{int a[N][N];node(){memset(a,0,sizeof a);}
}sqr;
node operator*(node x, node b){node c;for(int i=0;i<=n;i++){for(int j=0;j<=n;j++){for(int k=0;k<=n;k++)(c.a[i][j]+=x.a[i][k]*b.a[k][j]%mod)%=mod;}}return c;
}
void qpow(node &x, int y){node re;for(int i=0;i<=n;i++)re.a[i][i]=1;while(y){if(y&1)re=re*x;x=x*x;y>>=1;}x=re;
}
int main(){cin>>n>>m;for(int i=1,u,v;i<=m;i++)cin>>u>>v,sqr.a[u][v]=sqr.a[v][u]=1;cin>>t;for(int i=1;i<=n;i++)sqr.a[i][0]=1;for(int i=0;i<=n;i++)sqr.a[i][i]=1;qpow(sqr,t);int ans=0;for(int i=0;i<=n;i++)(ans+=sqr.a[1][i])%=mod;cout<<ans<<endl;}

P5343 【XR-1】分块

题面
方程很简单$dp[i]={\sum }_{j\in block,j\le i}dp[i-j]$，现在考虑如何矩阵优化。
由于块的大小不会超过$100$，所以我们开一个$100\times 100$的矩阵，首先预处理出$dp[1]-dp[100]$，将其填入$A$矩阵第一行中，再考虑所有$j\in block$，设$X[j][1]=1$，对于后99列设$X[j-1][j]=1$，则这样可以转移$A$矩阵，应用矩阵快速幂即可。

#include<bits/stdc++.h>
#define N 120
using namespace std;
const int mod=1e9+7;
const int len=100;
inline void read(int &x){int s=0,w=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}while(ch>='0'&&ch<='9'){s=(s<<3)+(s<<1)+(ch&15);ch=getchar();}x=s*w;
}
long long n;
int p,q,cnt,a[N],f[N],vis[N];
set<int> s;
struct node{int m[N][N];node(){memset(m,0,sizeof m);}
}sqr,A;
node operator*(node a, node b){node c;for(int i=0;i<=len;i++)for(int j=0;j<=len;j++)for(int k=0;k<=len;k++)c.m[i][j]=int((1LL*c.m[i][j]+1LL*a.m[i][k]*b.m[k][j]%mod)%mod);return c;
}
void qpow(node &x, long long y){node re;for(int i=0;i<=len;i++)re.m[i][i]=1;while(y){if(y&1)re=re*x;x=x*x;y>>=1;}x=re;
}
int main(){cin>>n;read(p);for(int i=1,x;i<=p;i++){read(x);if(s.find(x)==s.end())s.insert(x);}read(q);for(int i=1,x;i<=q;i++){read(x);if(s.find(x)!=s.end()&&!vis[x])a[++cnt]=x,vis[x]=1;}f[0]=1;for(int i=1;i<=len;i++)for(int j=1;j<=cnt;j++)if(a[j]<=i)f[i]=(1LL*f[i]+f[i-a[j]])%mod;if(n<=100){printf("%d\n",f[n]);return 0;}for(int i=0;i<=len;i++)sqr.m[0][len-i]=f[i];for(int i=1;i<=cnt;i++)A.m[a[i]-1][0]=1;for(int i=1;i<=len;i++)A.m[i-1][i]=1;qpow(A,n-100);sqr=sqr*A;printf("%d\n",sqr.m[0][0]);return 0;
}