这题怎么没人用矩乘啊
首先可以发现,这题的 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}
现在考虑倒着模拟亵渎的过程，可以发现相当于是要维护一个行向量集合，每次操作先给所有行向量乘上 $M$ ，然后加入一个初始行向量 $\{1,0,0...\}$ ，再给所有行向量乘 $M^b$ (其中 $b$ 这次亵渎删掉的段的右端点值)接着加入一个行向量 $\{-1,0,0...\}$，最后把所有行向量的第 $m+2$ 项加入答案。由于矩阵乘法具有分配律，所以可以直接维护所有行向量的和，然后进行一些矩阵加法和乘法即可做到 $O(Tm^4logn)$ 。注意到我们矩乘做的是行向量乘矩阵，所以单次矩乘可以 $O(m^2)$ ，然后我们再把转移矩阵的 $2^k$ 次幂预处理出来，就可以做到 $O(Tm^3logn)$。
大 力 出 奇 迹
```cpp
#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;
}
```