luogu P4593 [TJOI2018]教科书般的亵渎
这题怎么没人用矩乘啊
首先可以发现,这题的 \(k\) ,也就是亵渎使用次数为 \(m+1\) ，然后给出的 \(a_i\) 又会把 \([1,n]\) 划分成至多 \(m+1\) 个连续段。所以对每次亵渎，一个段 \([l,r]\) 给答案加上 \(\sum_{i=l}^r i^k\) ，并且所有连续段整体 \(-1\) ；此外如果某次亵渎后左端点最小的段的左端点为 \(0\) ，记它的右端点为 \(r_1\) ，那么会删掉这个段，并且其他段整体 \(-r_1\)。
那么我们模拟上述过程，然后每次对段操作前统计 \(\sum_{i=l}^r i^k\) 即可。考虑矩乘求这个东西，先把他拆成 \(\sum_{i=1}^r i^k-\sum_{i=1}^{l-1} i^k\) ，然后我们构造行向量 \(\vec{a}=\{x^0,x^1,x^2...x^{m+1},s\}\) ，即维护当前 \(i\) 的 \(0\) 到 \(m+1\) 次幂和答案，构造转移矩阵 \(M\) 就考虑 \((i+1)^k=\sum_{j=0}^{k}\binom{k}{j}i^j\) ，所以可以得到
- \(\forall 0\le i,j\le m+1,M_{i,j}=\binom{j}{i}\)
- \(M_{m+2,m+2}=1\)
- \(\forall 0\le i\le m+1,M_{i,m+2}=\binom{m+1}{i}\)
\(\vec{a}M^{n}\) 的第 \(m+2\) 项就是 \(\sum i^k\) 。不过注意到要做 \(m^2\) 次这样的求值，所以复杂度为 \(O(Tm^5logn)\)。
现在考虑倒着模拟亵渎的过程，可以发现相当于是要维护一个行向量集合，每次操作先给所有行向量乘上 \(M\) ，然后加入一个初始行向量 \(\{1,0,0...\}\) ，再给所有行向量乘 \(M^b\) (其中 \(b\) 这次亵渎删掉的段的右端点值)接着加入一个行向量 \(\{-1,0,0...\}\)，最后把所有行向量的第 \(m+2\) 项加入答案。由于矩阵乘法具有分配律，所以可以直接维护所有行向量的和，然后进行一些矩阵加法和乘法即可做到 \(O(Tm^4logn)\) 。注意到我们矩乘做的是行向量乘矩阵，所以单次矩乘可以 \(O(m^2)\) ，然后我们再把转移矩阵的 \(2^k\) 次幂预处理出来，就可以做到 \(O(Tm^3logn)\)。
大 力 出 奇 迹
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=53,mod=1e9+7;
const LL inf=8e18;
LL rd()
{
LL x=0,w=1;char ch=0;
while(ch<'0'||ch>'9'){if(ch=='-') w=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+(ch^48);ch=getchar();}
return x*w;
}
LL n,b[N],hp[N][2];
int m=51,t,p;
struct matrix
{
int a[N][N];
matrix(){memset(a,0,sizeof(a));}
void clr(){memset(a,0,sizeof(a));}
matrix operator + (const matrix &bb) const
{
matrix an;
for(int i=0;i<=m+1;++i)
for(int j=0;j<=m+1;++j)
an.a[i][j]=(a[i][j]+bb.a[i][j])%mod;
return an;
}
matrix operator * (const matrix &bb) const
{
matrix an;
for(int i=0;i<=m+1;++i)
for(int j=0;j<=m+1;++j)
{
LL nw=0;
for(int k=0;k<=m+1;++k)
{
nw+=1ll*a[i][k]*bb.a[k][j];
if(nw>inf) nw%=mod;
}
an.a[i][j]=nw%mod;
}
return an;
}
matrix operator & (const matrix &bb) const
{
matrix an;
for(int j=0;j<=m+1;++j)
{
LL nw=0;
for(int k=0;k<=m+1;++k)
{
nw+=1ll*a[0][k]*bb.a[k][j];
if(nw>inf) nw%=mod;
}
an.a[0][j]=nw%mod;
}
return an;
}
}aa,m1,bb[N][N];
int c[N][N];
void inii(int h)
{
for(int i=0;i<=h;++i)
for(int j=0;j<=h;++j)
bb[h][0].a[j][i]=c[i][j];
bb[h][0].a[h+1][h+1]=1;
for(int j=0;j<=h;++j)
bb[h][0].a[j][h+1]=c[h][j];
for(int i=1;i<N;++i) bb[h][i]=bb[h][i-1]*bb[h][i-1];
}
int main()
{
for(int i=0;i<N;++i)
{
c[i][0]=1;
for(int j=1;j<=i;++j) c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
}
aa.a[0][0]=1;
int T=rd();
while(T--)
{
n=rd(),m=rd();
b[0]=0;
for(int i=1;i<=m;++i) b[i]=rd();
sort(b+1,b+m+1),m=unique(b+1,b+m+1)-b;
while(m>1&&b[m-1]==n) --m,--n;
if(!bb[m][0].a[0][0]) inii(m);
b[m]=n+1,t=0;
for(int i=1;i<=m;++i)
{
if(b[i]==b[i-1]+1) continue;
hp[++t][0]=b[i-1]+1,hp[t][1]=b[i]-b[i-1]-1;
}
p=0;
while(t)
{
b[++p]=0;
for(int i=1;i<=t;++i) --hp[i][0];
if(!hp[1][0])
{
b[p]=hp[1][1];
for(int i=2;i<=t;++i) hp[i][0]-=hp[1][1];
for(int i=2;i<=t;++i) hp[i-1][0]=hp[i][0],hp[i-1][1]=hp[i][1];
--t;
}
}
m1.clr();
int ans=0;
while(p)
{
m1=m1*bb[m][0];
if(b[p])
{
m1.a[0][0]=(m1.a[0][0]+1)%mod;
for(int i=0;b[p];++i,b[p]>>=1)
if(b[p]&1) m1=m1&bb[m][i];
m1.a[0][0]=(m1.a[0][0]-1+mod)%mod;
}
ans=(ans+m1.a[0][m+1])%mod;
--p;
}
printf("%d\n",ans);
}
return 0;
}
