[loj6696]复读机 加强版
记$f(x)=\sum_{d\mid i}\frac{x^{i}}{i!}$，那么问题即求$n![x^{n}]f^{k}(x)$
记$\omega$为$d$次单位根，根据单位根反演有
$$
f(x)=\sum_{i\ge 0}\frac{\sum_{j=0}^{d-1}\omega^{ij}}{d}\frac{x^{i}}{i!}=\frac{1}{d}\sum_{j=0}^{d-1}\sum_{i\ge 0}\frac{(\omega^{j}x)^{i}}{i!}=\frac{1}{d}\sum_{j=0}^{d-1}e^{\omega^{j}x}
$$
根据$\omega^{2}=-1$（$d=4$时）和$\omega^{2}=\omega-1$（$d=6$时）降幂，将$k$次方展开后，每一项即形如$e^{(a+b\omega )x}$
（提出$\frac{1}{d^{k}}$的常数）记$S_{a,b}$为$e^{(a +b\omega)x}$的系数，问题也即求$\frac{1}{d^{k}}\sum_{|a|,|b|\le k}(a+b\omega)^{n}S_{a,b}$，快速幂计算$(a+b\omega)^{n}$即可做到$o(k^{2}\log n)$
另外，这里的$\omega$并不需要扩域，注意到$19491001$有原根$g=7$，进而模意义下$\omega=g^{\frac{mod-1}{d}}$
关于$S_{a,b}$，设$\omega^{j}=a_{j}+b_{j}\omega,G(x,y)=\sum_{j=0}^{d-1}x^{a_{j}}y^{b_{j}}$，则$S$的二元生成函数为$S(x,y)=G^{k}(x,y)$
将其简单变形，不难得到$[x^{a}y^{b}]G(x,y)S'(x,y)=k\cdot [x^{a}y^{b}]G'(x,y)S(x,y)$（求导均指对$x$的偏导）
当$d=4$时，简单计算可得$\begin{cases}G(x,y)=x+y+x^{-1}+y^{-1}\\G'(x,y)=1-x^{-2}\end{cases}$，代入可得
$$
aS_{a,b}+(a+1)S_{a+1,b-1}+(a+2)S_{a+2,b}+(a+1)S_{a+1,b+1}=k(S_{a,b}-S_{a+2,b})
$$
当$d=6$时，简单计算可得$\begin{cases}G(x,y)=x+y+x^{-1}y+x^{-1}+y^{-1}+xy^{-1}\\G'(x,y)=1-x^{-2}y-x^{-2}+y^{-1}\end{cases}$，同理代入可得
$$
aS_{a,b}+(a+1)S_{a+1,b-1}+(a+2)S_{a+2,b-1}+(a+2)S_{a+2,b}+(a+1)S_{a+1,b+1}+aS_{a,b+1}=k(S_{a,b}-S_{a+2,b-1}-S_{a+2,b}+S_{a,b+1})
$$
两者均可$o(k^{2})$线性递推（初始状态为$S_{-k,*}$，注意乘$k!$），总复杂度为$o(k^{2}\log n)$，可以通过
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define N 4005 4 #define mod 19491001 5 #define g 7 6 #define ll long long 7 int n,k,d,w,ans,inv[N],Inv[N],S[N][N]; 8 int get_inv(int x){ 9 if (x>0)return inv[x]; 10 return mod-inv[-x]; 11 } 12 int C(int n,int m){ 13 int s=1; 14 for(int i=1;i<=m;i++)s=(ll)s*(n-i+1)%mod*inv[i]%mod; 15 return s; 16 } 17 int qpow(int n,int m){ 18 int s=n,ans=1; 19 while (m){ 20 if (m&1)ans=(ll)ans*s%mod; 21 s=(ll)s*s%mod,m>>=1; 22 } 23 return ans; 24 } 25 int main(){ 26 inv[0]=inv[1]=1; 27 for(int i=2;i<N;i++)inv[i]=(ll)(mod-mod/i)*inv[mod%i]%mod; 28 scanf("%d%d%d",&n,&k,&d),w=qpow(g,(mod-1)/d); 29 if (d==4)S[0][k]=1; 30 else{ 31 for(int j=0;j<=k;j++)S[0][j+k]=C(k,j); 32 } 33 for(int i=-k+1;i<=k;i++) 34 for(int j=-k;j<=k;j++){ 35 int x=i+k,y=j+k; 36 if (d==4){ 37 if (x>1)S[x][y]=(S[x][y]+(ll)(k-i+2)*S[x-2][y])%mod; 38 if ((x)&&(y))S[x][y]=(S[x][y]+(ll)(mod-i+1)*S[x-1][y-1])%mod; 39 if (x)S[x][y]=(S[x][y]+(ll)(mod-i+1)*S[x-1][y+1])%mod; 40 S[x][y]=(ll)get_inv(k+i)*S[x][y]%mod; 41 } 42 if (d==6){ 43 if (x>1)S[x][y]=(S[x][y]+(ll)(k-i+2)*S[x-2][y])%mod; 44 if ((x)&&(y))S[x][y]=(S[x][y]+(ll)(mod-i+1)*S[x-1][y-1])%mod; 45 if (y)S[x][y]=(S[x][y]-(ll)(k+i)*S[x][y-1]%mod+mod)%mod; 46 if (x)S[x][y]=(S[x][y]+(ll)(mod-i+1)*S[x-1][y+1])%mod; 47 if (x>1)S[x][y]=(S[x][y]+(ll)(k-i+2)*S[x-2][y+1])%mod; 48 S[x][y]=(ll)get_inv(k+i)*S[x][y]%mod; 49 } 50 } 51 for(int i=-k;i<=k;i++) 52 for(int j=-k;j<=k;j++) 53 if (S[i+k][j+k])ans=(ans+(ll)qpow((i+(ll)j*w)%mod+mod,n)*S[i+k][j+k])%mod; 54 ans=(ll)ans*qpow(get_inv(d),k)%mod; 55 printf("%d\n",ans); 56 return 0; 57 }
