这题怎么没人用矩乘啊
首先可以发现,这题的 (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}inom{k}{j}i^j) ，所以可以得到
- (forall 0le i,jle m+1,M_{i,j}=inom{j}{i})
- (M_{m+2,m+2}=1)
- (forall 0le ile m+1,M_{i,m+2}=inom{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
",ans);
}
return 0;
}