应队友要求,开始学线性代数,具体路线是矩阵 → \rightarrow →高斯消元 → \rightarrow →线性基。为多项式做个准备
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 = [ a 3 a 2 a 1 0 0 0 0 0 0 ] , X = [ 1 1 0 0 0 1 1 0 0 ] , A_3=\left [ \begin{matrix} a_3& a_2 & a_1 \\ 0& 0 &0 \\ 0 & 0 & 0 \\ \end{matrix} \right] ,X=\left [ \begin{matrix} 1& 1 & 0 \\ 0& 0 &1 \\ 1& 0 & 0 \\ \end{matrix} \right], A3= a300a200a100 ,X= 101100010 ,
则 A 3 X = [ a 4 a 3 a 2 0 0 0 0 0 0 ] = A 4 , A_3X=\left [ \begin{matrix} a_4& a_3 & a_2 \\ 0& 0 &0 \\ 0 & 0 & 0 \\ \end{matrix} \right]=A_4, A3X= a400a300a200 =A4,
因此对 X X 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 = [ f 2 f 1 0 0 ] , X = [ 1 1 1 0 ] , A_2=\left [ \begin{matrix} f_2 & f_1 \\ 0 &0 \\ \end{matrix} \right] ,X=\left [ \begin{matrix} 1& 1 \\ 1& 0 \\ \end{matrix} \right], A2=[f20f10],X=[1110],
则 A 2 X = [ f 3 f 2 0 0 ] = A 3 , A_2X=\left [ \begin{matrix} f_3 & f_2 \\ 0 &0 \\ \end{matrix} \right]=A_3, A2X=[f30f20]=A3,
#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 = [ f 2 f 1 0 0 ] , X = [ P 1 Q 0 ] , A_2=\left [ \begin{matrix} f_2 & f_1 \\ 0 &0 \\ \end{matrix} \right] ,X=\left [ \begin{matrix} P& 1 \\ Q& 0 \\ \end{matrix} \right], A2=[f20f10],X=[PQ10],
则 A 2 X = [ f 3 f 2 0 0 ] = A 3 , A_2X=\left [ \begin{matrix} f_3 & f_2 \\ 0 &0 \\ \end{matrix} \right]=A_3, A2X=[f30f20]=A3,
#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 18 18年外出集训堵了个论文费了两三天时间才切了的人生中第一道黑题,现在变成紫题了。
有一个性质是 f n m o d p f_n\mod p fnmodp有循环节,且循环节长度不会超过 6 p 6p 6p,这还有个名叫皮萨诺定理。所以我们考虑求出循环节的长度,然后用矩阵乘法求出结果。
引理:对于 f n m o d p f_n\mod p fnmodp的循环节 g ( p ) g(p) g(p)有如下性质:
1. p = p i α i 1.p=p_i^{\alpha_i} 1.p=piαi,即 p p p为质数的幂时, g ( p ) = g ( p i ) × p i α i − 1 g(p)=g(p_i)\times p_i^{\alpha_i-1} g(p)=g(pi)×piαi−1
2. p = ∏ p i α i 2.p=\prod p_i^{\alpha_i} 2.p=∏piαi,即 p p p为合数时, g ( p ) = l c m ( g ( p i α i ) g(p)=lcm(g(p_i^{\alpha_i}) g(p)=lcm(g(piαi)
对于 g ( p ) g(p) g(p)这么算,如果 5 5 5是模 p p p的二次剩余,那么循环节为 p − 1 p-1 p−1的因子,否则为 2 p + 2 2p+2 2p+2的因子。
因为 p p p不是特别大,直接取 p − 1 p-1 p−1和 2 p + 2 2p+2 2p+2即可。
对于 p ≤ 5 p\le 5 p≤5就暴力算即可, g ( 2 ) = 3 , g ( 3 ) = 5 , g ( 5 ) = 20 g(2)=3,g(3)=5,g(5)=20 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】分块
题面
方程很简单 d p [ i ] = ∑ j ∈ b l o c k , j ≤ i d p [ i − j ] dp[i]=\sum_{j\in block,j\le i}dp[i-j] dp[i]=∑j∈block,j≤idp[i−j],现在考虑如何矩阵优化。
由于块的大小不会超过 100 100 100,所以我们开一个 100 × 100 100\times 100 100×100的矩阵,首先预处理出 d p [ 1 ] − d p [ 100 ] dp[1]-dp[100] dp[1]−dp[100],将其填入 A A A矩阵第一行中,再考虑所有 j ∈ b l o c k j\in block j∈block,设 X [ j ] [ 1 ] = 1 X[j][1]=1 X[j][1]=1,对于后99列设 X [ j − 1 ] [ j ] = 1 X[j-1][j]=1 X[j−1][j]=1,则这样可以转移 A A 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;
}