NFLSOJ #10319. -「2020联考北附2」碎梦（NTT+对角化）
作者：互联网
一个奇怪的解法，时间复杂度比标算少一个 \(k^2m\) 所以建议把数据范围加强到 \(k=300,m=3\times 10^4\)，8s 时限把标算卡掉（bushi）
首先考虑部分分解法。对于 \(m=100\)，注意到转移可以写成矩阵的形式，具体来说记 \(A=\begin{bmatrix}0&1&0&0&\cdots&0&0\\0&1&1&0&\cdots&0&0\\0&0&2&1&\cdots&0&0\\0&0&0&3&\cdots&0&0\\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&0&\cdots&m-1&1\\0&0&0&0&\cdots&0&m\end{bmatrix}\)，那么有 \(\begin{bmatrix}f_{i,0}&f_{i,1}&\cdots&f_{i,m}\end{bmatrix}=\begin{bmatrix}f_{i-1,0}&f_{i-1,1}&\cdots&f_{i-1,m}\end{bmatrix}\times A\)。因此我们可以仿照 CF576D 或者 P6772 的套路，将全部 \(k\) 个三元组排个序（事实上输入已经有序了，所以不必再 sort），然后每次用一个 \(1\times(m+1)\) 的矩阵维护当前的 \(f_i\)，对于一个特殊点 \((x_i,y_i,z_i)\)，就矩阵快速幂推到对应的 \(x_i\) 处然后将 \(y_i\) 位置上的值改为 \(z_i\) 即可。时间复杂度 \(km^3\log n\)。
考虑优化。注意到该矩阵事上三角矩阵且对角线上元素刚好是 \(0,1,2,\cdots,m\)，这让我们不禁想起了上周二模拟赛的 T4，因此考虑按照那题的套路对角化。具体来说我们 \(m^2\) 地求出矩阵的特征向量然后拼起来得到矩阵 \(T\)，再一遍求逆求出 \(T\) 的逆矩阵，这样根据对角化那一套理论 \(A^p=TB^pT^{-1}\)，其中 \(B_{i,i}=i,B_{i,j}=0(i\ne j)\)。由于矩阵乘向量复杂度可实现 \(\mathcal O(m^2)\)。复杂度可实现 \(km^2+m^3+km\log n\)，可以获得 85 分的好成绩。
继续优化，打个表发现矩阵 \(T\) 很特别，具体来说，
\[T_{i,j}=\begin{cases}0&(i<j)\\\dfrac{1}{(i-j)!}&(i\ge j)\end{cases} \]也就是说对于一个 \(1\times (m+1)\) 的矩阵 \(X\)，设 \(XT=Y\)，那么由于所有 \(i-j\) 相同的 \(T_{i,j}\) 都相同，因此转移可以写成卷积的形式，即 \(Y_i=\sum\limits_{j=0}^iX_j\dfrac{1}{(i-j)!}\)，因此一个向量乘 \(T\) 可用 NTT 优化至 \(m\log m\)。同理 \(T^{-1}\) 也有类似的性质，即
\[T^{-1}_{i,j}=\begin{cases}0&(i<j)\\\dfrac{1}{(i-j)!}(-1)^{i-j}&(i\ge j)\end{cases} \]也可 NTT 优化。而 \(B^p\) 是对角矩阵，乘法自然也可 \(\mathcal O(m)\) 算出，再加上快速幂的 \(\log\)，总复杂度也是 \(m\log n\)，这样向量乘 \(A^p\) 就被成功优化至 \(m\log n\)，总复杂度 \(km\log n\)，可以通过。
关于卡常：现场写的这个复杂度的做法被卡成了 90 分，原因是在向量乘 \(A^p\) 的过程中要 6 遍 NTT，不过注意到有两次 NTT 过程可以预处理出来，这样每次只要 4 次 NTT，还有就是 NTT 本身板子的优化，感谢 chasedeath 神仙提供常数更小的 NTT 板子 %%%%%%%
int n,m,k;
int qpow(int x,int e){
int ret=1;
for(;e;e>>=1,x=1ll*x*x%MOD) if(e&1) ret=1ll*ret*x%MOD;
return ret;
}
namespace sub1{
const int MAXM=500;
int a[MAXM+5][MAXM+5],b[MAXM+5][MAXM+5],c[MAXM+5][MAXM+5],d[MAXM+5][MAXM*2+5];
void calcc(int x){
for(int i=0;i<=m;i++) for(int j=0;j<=m;j++) b[i][j]=a[i][j];
for(int i=0;i<=m;i++) b[i][i]=(b[i][i]-x+MOD)%MOD;
c[x][x]=1;
for(int i=x-1;~i;i--){
c[i][x]=0;
for(int j=i+1;j<=x;j++) c[i][x]=(c[i][x]-1ll*c[j][x]*b[i][j]%MOD+MOD)%MOD;
c[i][x]=1ll*c[i][x]*qpow(b[i][i],MOD-2)%MOD;
}
}
void calcd(){
for(int i=0;i<=m;i++) for(int j=0;j<=m;j++) d[i][j]=c[i][j];
for(int i=0;i<=m;i++) d[i][i+m+1]=1;
for(int i=m;~i;i--){
for(int j=i+1;j<=m;j++){
for(int k=m+1;k<=2*m+1;k++)
d[i][k]=(d[i][k]-1ll*d[i][j]*d[j][k]%MOD+MOD)%MOD;
}
}
for(int i=0;i<=m;i++) for(int j=0;j<=m;j++) swap(d[i][j],d[i][j+m+1]);
}
struct mat{
int n,m,a[MAXM+5][MAXM+5];
mat(){n=m=0;memset(a,0,sizeof(a));}
mat operator *(const mat &rhs){
mat ret;ret.n=n;ret.m=rhs.m;
for(int i=0;i<=n;i++) for(int j=0;j<=rhs.m;j++)
for(int k=0;k<=m;k++) ret.a[i][j]=(ret.a[i][j]+1ll*a[i][k]*rhs.a[k][j])%MOD;
return ret;
}
} A,_A,cur,trs;
int solve(){
for(int i=0;i<=m;i++) a[i][i]=i;
for(int i=1;i<=m;i++) a[i-1][i]=1;
for(int i=0;i<=m;i++) calcc(i);
calcd();A.n=A.m=_A.n=_A.m=m;
for(int i=0;i<=m;i++) for(int j=0;j<=m;j++){
A.a[i][j]=c[i][j];
_A.a[i][j]=d[i][j];
}
// for(int i=0;i<=m;i++) for(int j=0;j<=m;j++)
// printf("%d%c",_A.a[i][j]," \n"[j==m]);
cur.n=0;cur.m=m;cur.a[0][0]=1;
int pre=0;
for(int i=1;i<=k;i++){
int x,y,z;scanf("%d%d%d",&x,&y,&z);
int len=x-pre;
cur=cur*A;
for(int j=0;j<=m;j++) cur.a[0][j]=1ll*cur.a[0][j]*qpow(j,len)%MOD;
cur=cur*_A;pre=x;
cur.a[0][y]=z;
} int len=n-pre;cur=cur*A;
for(int j=0;j<=m;j++) cur.a[0][j]=1ll*cur.a[0][j]*qpow(j,len)%MOD;
cur=cur*_A;printf("%d\n",cur.a[0][m]);
return 0;
}
}
namespace sub2{
const int MAXP=1<<18;
const int pr=3;
const int ipr=332748118;
int rev[MAXP+5];
void NTT(int *a,int len,int type){
int lg=31-__builtin_clz(len);
for(int i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<lg-1);
for(int i=0;i<len;i++) if(rev[i]<i) swap(a[i],a[rev[i]]);
static int e[MAXP+5];e[0]=1;
for(int i=2;i<=len;i<<=1){
int W=qpow((type<0)?ipr:pr,(MOD-1)/i);
for(int j=(i>>1)-2;j>=0;j-=2) e[j+1]=1ll*(e[j]=e[j>>1])*W%MOD;
for(int j=0;j<len;j+=i){
for(int k=0;k<(i>>1);k++){
int X=a[j+k],Y=1ll*a[(i>>1)+j+k]*e[k]%MOD;
a[j+k]=X+Y>=MOD?X+Y-MOD:X+Y;
a[(i>>1)+j+k]=X-Y<0?X-Y+MOD:X-Y;
}
}
}
if(!~type){
int iv=qpow(len,MOD-2);
for(int i=0;i<len;i++) a[i]=1ll*a[i]*iv%MOD;
}
}
int A[MAXP+5],B[MAXP+5],LEN=1;
int cur[MAXP+5],ifac[MAXP+5],ifac_[MAXP+5],ifac__[MAXP+5];
void times(int x){
memset(A,0,sizeof(A));
for(int i=0;i<=m;i++) A[i]=cur[i];
NTT(A,LEN,1);
for(int i=0;i<LEN;i++) A[i]=1ll*A[i]*ifac_[i]%MOD;
NTT(A,LEN,-1);
for(int i=0;i<=m;i++) cur[i]=1ll*A[i]*qpow(i,x)%MOD;
memset(A,0,sizeof(A));
for(int i=0;i<=m;i++) A[i]=cur[i];
NTT(A,LEN,1);
for(int i=0;i<LEN;i++) A[i]=1ll*A[i]*ifac__[i]%MOD;
NTT(A,LEN,-1);
for(int i=0;i<=m;i++) cur[i]=A[i];
}
int solve(){
for(int i=(ifac[0]=ifac[1]=1)+1;i<=MAXP;i++) ifac[i]=1ll*ifac[MOD%i]*(MOD-MOD/i)%MOD;
for(int i=1;i<=MAXP;i++) ifac[i]=1ll*ifac[i-1]*ifac[i]%MOD;
while(LEN<=m+m) LEN<<=1;cur[0]=1;int pre=0;
for(int i=0;i<=m;i++) ifac_[i]=ifac[i];NTT(ifac_,LEN,1);
for(int i=0;i<=m;i++) ifac__[i]=(i&1)?(MOD-ifac[i]):ifac[i];
NTT(ifac__,LEN,1);
for(int i=1;i<=k;i++){
int x,y,z;scanf("%d%d%d",&x,&y,&z);
times(x-pre);pre=x;cur[y]=z;
} times(n-pre);printf("%d\n",cur[m]);
return 0;
}
}
int main(){
freopen("broken.in","r",stdin);
freopen("broken.out","w",stdout);
scanf("%d%d%d",&n,&m,&k);
if(m<=100) return sub1::solve();
return sub2::solve();
}
/*
900000000 100000 20
5 2 1
10 5 2
50 20 10
100 50 20
500 200 100
1000 500 200
5000 2000 1000
10000 5000 2000
50000 20000 10000
100000 50000 10000
500000 100000 100000
1000000 100000 1000000
5000000 100000 10000000
50000000 100000 100000000
100000000 100000 200000000
200000000 100000 300000000
300000000 100000 400000000
400000000 100000 500000000
500000000 100000 600000000
600000000 100000 700000000
*/
标签：10319,北附,int,复杂度,矩阵,NTT,cdots,MAXM,联考 来源： https://www.cnblogs.com/ET2006/p/NFLSOJ-10319.html