线性代数都快忘完了,先挖个坑...
矩阵加速递推
我们将第n项斐波那契数列记作
考虑 ,把斐波那契数列的相邻两项放在一个行矩阵里面,就像这个样子:
我们需要把 变成
需要用到矩阵乘法
因为我们一开始的是
所以我们只需要将其乘以
即可得到
当n很大的时候比如为1e18,通过直接乘是肯定会TLE的,所以考虑使用快速幂,这样就可以在log的时间里面求出答案,加速成功
#include <cstdio>
using namespace std;
struct Matrix{
static const int MOD = 10007;
int data[2][2];
Matrix(int data00, int data01, int data10, int data11){
data[0][0] = data00; data[0][1] = data01;
data[1][0] = data10; data[1][1] = data11;
}
Matrix operator * (const Matrix &a) const{
Matrix ans = Matrix(0, 0, 0, 0);
for(int i = 0; i<2; ++i)
for(int j = 0; j<2; ++j)
for(int k = 0; k<2; ++k)
ans.data[i][j] = (ans.data[i][j] + data[i][k]*a.data[k][j])%MOD;
return ans;
}
Matrix pow( int power){
Matrix ans = Matrix(1, 0, 0, 1);
Matrix base = (*this);
while(power){
if(power&1) ans = ans * base;
base = base * base;
power >>= 1;
}
return ans;
}
};
int main(){
int n;
scanf("%d", &n);
Matrix ans = Matrix(1, 1, 1, 0);
printf("%d\n", ans.pow(n).data[1][0]);
return 0;
}
线段树维护矩阵乘法
把1-6的操作用乘上对应的矩阵实现,那么只要使用线段树维护矩阵即可
#include<bits/stdc++.h>
using namespace std;
#define int long long
const int maxn=2.5e5+5;
const int mod=998244353;
struct Matrix{
int n,m;
int a[5][5];
Matrix(){}
Matrix(int _n,int _m):n(_n),m(_m){
memset(a,0,sizeof(a));
}
Matrix operator+(const Matrix& rhs){
Matrix res(n,m);
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
res.a[i][j]=(a[i][j]+rhs.a[i][j])%mod;
}
}
return res;
}
Matrix operator*(const Matrix &rhs){
Matrix res(n,rhs.m);
for(int i=1;i<=n;i++){
for(int j=1;j<=m;j++){
for(int k=1;k<=rhs.m;k++){
res.a[i][k]=(res.a[i][k]+a[i][j]*rhs.a[j][k])%mod;
}
}
}
return res;
}
};
Matrix num[maxn];
struct SegTree{
int l,r;
Matrix sum,tag;
}tr[maxn<<2];
Matrix base,base1,base2,base3,base4,base5,base6;
void init(){
base=Matrix(4,4);
for(int i=1;i<=4;i++){
base.a[i][i]=1;
}
/*
base1:| base2:| base3:
1 0 0 0| 1 0 0 0| 1 0 1 0
1 1 0 0| 0 1 0 0| 0 1 0 0
0 0 1 0| 0 1 1 0| 0 0 1 0
0 0 0 1| 0 0 0 1| 0 0 0 1
*/
base1=base,base2=base,base3=base;
base1.a[2][1]=1;
base2.a[3][2]=1;
base3.a[1][3]=1;
/*
base4: | base5: | base6:
1 0 0 0 | 1 0 0 0 | 1 0 0 0
0 1 0 0 | 0 k 0 0 | 0 1 0 0
0 0 1 0 | 0 0 1 0 | 0 0 0 0
k 0 0 1 | 0 0 0 1 | 0 0 k 1
*/
base4=base,base5=base,base6=base;
base6.a[3][3]=0;
}
void build(int p,int l,int r){
tr[p].l=l,tr[p].r=r;
if(l==r){
tr[p].sum=num[l];
tr[p].tag=base;
return;
}
int m=(l+r)>>1;
build(p<<1,l,m);
build(p<<1|1,m+1,r);
tr[p].sum=tr[p<<1].sum+tr[p<<1|1].sum;
tr[p].tag=base;
return;
}
void pushdown(int p){
auto &root=tr[p],&left=tr[p<<1],&right=tr[p<<1|1];
left.sum=left.sum*root.tag;
right.sum=right.sum*root.tag;
left.tag=left.tag*root.tag;
right.tag=right.tag*root.tag;
root.tag=base;
return;
}
void segadd(int p,int l,int r,Matrix k){
if(tr[p].l>=l&&tr[p].r<=r){
tr[p].tag=tr[p].tag*k;
tr[p].sum=tr[p].sum*k;
return;
}
pushdown(p);
int m=(tr[p].l+tr[p].r)>>1;
if(l<=m){
segadd(p<<1,l,r,k);
}
if(r>m){
segadd(p<<1|1,l,r,k);
}
tr[p].sum=tr[p<<1].sum+tr[p<<1|1].sum;
return;
}
Matrix segsum(int p,int l,int r){
if(tr[p].l>=l&&tr[p].r<=r){
return tr[p].sum;
}
pushdown(p);
int m=(tr[p].l+tr[p].r)>>1;
Matrix res(4,4);
if(l<=m){
res=res+segsum(p<<1,l,r);
}
if(r>m){
res=res+segsum(p<<1|1,l,r);
}
return res;
}
signed main(){
init();
int n;
cin>>n;
for(int i=1;i<=n;i++){
num[i]=Matrix(1,4);
for(int j=1;j<=3;j++){
cin>>num[i].a[1][j];
}
num[i].a[1][4]=1;
}
build(1,1,n);
int m;
cin>>m;
while(m--){
int op,l,r,k;
cin>>op>>l>>r;
if(op==1){
segadd(1,l,r,base1);
}
else if(op==2){
segadd(1,l,r,base2);
}
else if(op==3){
segadd(1,l,r,base3);
}
else if(op==4){
cin>>k;
base4.a[4][1]=k;
segadd(1,l,r,base4);
}
else if(op==5){
cin>>k;
base5.a[2][2]=k;
segadd(1,l,r,base5);
}
else if(op==6){
cin>>k;
base6.a[4][3]=k;
segadd(1,l,r,base6);
}
else if(op==7){
Matrix res=segsum(1,l,r);
cout<<res.a[1][1]<<" "<<res.a[1][2]<<" "<<res.a[1][3]<<"\n";
}
}
}
高斯消元 (求行列式,时间复杂度n^3)
const double EPS = 1E-9;
int n;
vector<vector<double> > a(n, vector<double>(n));
double det = 1;
for (int i = 0; i < n; ++i) {
int k = i;
for (int j = i + 1; j < n; ++j)
if (abs(a[j][i]) > abs(a[k][i])) k = j;
if (abs(a[k][i]) < EPS) {
det = 0;
break;
}
swap(a[i], a[k]);
if (i != k) det = -det;
det *= a[i][i];
for (int j = i + 1; j < n; ++j) a[i][j] /= a[i][i];
for (int j = 0; j < n; ++j)
if (j != i && abs(a[j][i]) > EPS)
for (int k = i + 1; k < n; ++k) a[j][k] -= a[i][k] * a[j][i];
}
cout << det;
双指针规避矩阵除法
对于一个矩阵连续区间l-r相乘,如果要求l+1-r,即删掉l这个矩阵,那么就要乘以l这个矩阵的逆矩阵,这个可以用高斯消元解决,但是如果矩阵不是可逆的,那么就不能执行删除操作了,所以考虑使用双指针的办法l,r用来表示当前的左右端点,如果可行则r++,否则l++,但是这样还是没有实现删除操作,所以另外加一个limit,base矩阵代表limit+1-r的乘积,b[i]预处理一下,代表i-limit的乘积,当目前可行时,则r++,否则l++,如果l>limit,那么就把limit赋值为r,然后预处理算出limit+1-r的b[i],并且令因为当前limit+1=r+1,limit+1-r目前没有矩阵,所以把base赋值为1,重新开始新的一轮,这样可以实现O(1)的转移
#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int inf=1e9+8;
ll n,m,k;
struct square{
ll v[22][22];
square(int x=0){
memset(v,0,sizeof(v));
for(int i=1;i<=m;i++){
v[i][i]=x;
}
}
square operator *(square q){
square ans;
for(int i=1;i<=m;i++){
for(int j=1;j<=m;j++){
for(int k=1;k<=m;k++){
ans.v[i][j]=ans.v[i][j]+v[i][k]*q.v[k][j];
if(ans.v[i][j]>inf){//多次累乘可能越界,过大就处理一下
ans.v[i][j]=inf;
}
}
}
}
return ans;
}
void print(){
for(int i=1;i<=m;i++){
for(int j=1;j<=m;j++){
cout<<v[i][j]<<" ";
}
cout<<endl;
}
}
};
square a[5010];
square b[5010];
bool check(square x,square y){//判断当前方案数是否超出k
ll ans=0;
for(int i=1;i<=m;i++){
ans+=x.v[1][i]*y.v[i][m];
if(ans>k){
return 0;
}
}
return 1;
}
int main(){
ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);//同步流
int t;
cin>>t;
while(t--){
cin>>n>>m>>k;
for(int i=1;i<=n;i++){
int l;
cin>>l;
a[i]=1;
for(int j=1;j<=l;j++){
int u,v;
cin>>u>>v;
a[i].v[u][v]=1;
}
}
int ans=0;
b[0].v[1][m]=inf;//初始化
square base=1;
for(int l=0,lim=0,r=1;r<=n;r++){
base=base*a[r];//储存a[lim+1]-a[r]的累乘和
while(!check(b[l],base)){//方案数过多,开始缩小区间,查找当前人对应的l
//根据base,和b[i]的含义,b[l]*base即为a[l]-a[r]的累乘和
l++;
if(l>lim){//l超过lim了,此时lim后面的b还没有定义,需要更新lim和对应的b
base=1;
b[r]=a[r];//b[i]储存的是a[i]-a[lim]间的累乘和,
for(int i=r-1;i>lim;i--){
b[i]=a[i]*b[i+1];
}
lim=r;
}
}
ans=max(ans,r-l+1);//更新答案
}
cout<<ans<<endl;
}
return 0;
}