[集训队作业2018] 复读机（生成函数，单位根反演）

传送门

subtask 1：d=1d=1d=1

答案为knk^nkn。

subtask 2：n≤1000,k≤100n\leq1000,k\leq 100n≤1000,k≤100

设f[i][j]f[i][j]f[i][j]表示由iii个复读机来分jjj个时间点的方案数。
可以得到递推式：
f[i][j]=∑p=0j[d∣p]Cjp×f[i−1][j−p]f[i][j]=\sum_{p=0}^{j}[d|p]C_{j}^{p}\times f[i-1][j-p]f[i][j]=p=0∑j[d∣p]Cjp×f[i−1][j−p]
O(nk2)O(nk^2)O(nk2)暴力DP即可。

subtask 3：d=2d=2d=2

把d=2d=2d=2代入上面的递推式得：
f[i][j]=∑p=0j[2∣p]Cjp×f[i−1][j−p]f[i][j]=\sum_{p=0}^{j}[2|p]C_{j}^{p}\times f[i-1][j-p]f[i][j]=p=0∑j[2∣p]Cjp×f[i−1][j−p]
f[i][j]=∑p=0j[2∣p]j!p!(j−p)!×f[i−1][j−p]f[i][j]=\sum_{p=0}^{j}[2|p]\frac{j!}{p!(j-p)!}\times f[i-1][j-p]f[i][j]=p=0∑j[2∣p]p!(j−p)!j!×f[i−1][j−p]
f[i][j]j!=∑p=0j[2∣p]1p!f[i−1][j−p](j−p)!\frac{f[i][j]}{j!}=\sum_{p=0}^{j}[2|p]\frac{1}{p!}\frac{f[i-1][j-p]}{(j-p)!}j!f[i][j]=p=0∑j[2∣p]p!1(j−p)!f[i−1][j−p]
设Ai(x)=∑j=0∞f[i][j]×xjj!,B(x)=∑j=0∞[2∣j]xjj!A_i(x)=\sum_{j=0}^{\infty}f[i][j]\times\frac{x^j}{j!},B(x)=\sum_{j=0}^{\infty}[2|j]\frac{x^j}{j!}Ai(x)=∑j=0∞f[i][j]×j!xj,B(x)=∑j=0∞[2∣j]j!xj
代入上式得：
Ai(x)[xj]=∑p=0jB(x)[xp]×Ai−1(x)[xj−p]A_i(x)[x^j]=\sum_{p=0}^{j}B(x)[x^p]\times A_{i-1}(x)[x^{j-p}]Ai(x)[xj]=p=0∑jB(x)[xp]×Ai−1(x)[xj−p]
所以Ai(x)=B(x)Ai−1(x)A_i(x)=B(x)A_{i-1}(x)Ai(x)=B(x)Ai−1(x)
又因为A0(x)=1A_0(x)=1A0(x)=1
所以Ai(x)=Bi(x)A_i(x)=B^i(x)Ai(x)=Bi(x)

化简B(x)B(x)B(x)：
因为ex=∑j=0∞xjj!,e−x=∑j=0∞(−1)j×xjj!e^x=\sum_{j=0}^{\infty}\frac{x^j}{j!},e^{-x}=\sum_{j=0}^{\infty}(-1)^j\times \frac{x^j}{j!}ex=∑j=0∞j!xj,e−x=∑j=0∞(−1)j×j!xj
所以B(x)=ex+e−x2B(x)=\frac{e^x+e^{-x}}{2}B(x)=2ex+e−x

所以f[k][n]f[k][n]f[k][n](最终答案)为(ex+e−x2)k×n!(\frac{e^x+e^{-x}}{2})^k\times n!(2ex+e−x)k×n!的nnn次项系数。
将上式二项式展开得：n!×12k∑i=0kCki×eix×e(i−k)xn!\times\frac{1}{2^k}\sum_{i=0}^{k}C_{k}^{i}\times e^{ix}\times e^{(i-k)x}n!×2k1i=0∑kCki×eix×e(i−k)x
=n!×12k∑i=0kCki×e(2i−k)x=n!\times\frac{1}{2^k}\sum_{i=0}^{k}C_{k}^{i}\times e^{(2i-k)x}=n!×2k1i=0∑kCki×e(2i−k)x
而e(2i−k)x=∑j=0∞((2i−k)x)jj!e^{(2i-k)x}=\sum_{j=0}^{\infty}\frac{((2i-k)x)^{j}}{j!}e(2i−k)x=∑j=0∞j!((2i−k)x)j，其nnn次项系数为(2i−k)nn!\frac{(2i-k)^n}{n!}n!(2i−k)n

所以最终答案为12k∑i=0kCki×(2i−k)n\frac{1}{2^k}\sum_{i=0}^{k}C_{k}^{i}\times (2i-k)^n2k1i=0∑kCki×(2i−k)n

subtask 4：d=3

同样考虑生成函数，答案就是(∑i=0∞[3∣i]xii!)k×n!(\sum_{i=0}^{\infty}[3|i]\frac{x^i}{i!})^k\times n!(∑i=0∞[3∣i]i!xi)k×n!的nnn次项系数。

有一个trick，叫单位根反演，大概是这样：
[n∣k]=1n∑i=0n−1wnki[n|k]=\frac{1}{n}\sum_{i=0}^{n-1}w_{n}^{ki}[n∣k]=n1i=0∑n−1wnki
注意到19491001−119491001-119491001−1是3的倍数，即mod−1mod-1mod−1是3的倍数，故存在三次单位根rrr。由单位根的定义知，r=Rmod−13r=R^{\frac{mod-1}{3}}r=R3mod−1 ，其中RRR是modmodmod的一个原根。运用单位根反演，有：

∑i=0∞[3∣i]xii!=13∑i=0∞(r0+ri+r2i)xii!=ex+erx+er2x3\sum_{i=0}^{\infty}[3|i]\frac{x^i}{i!}=\frac{1}{3}\sum_{i=0}^{\infty}\frac{(r^0+r^{i}+r^{2i})x^i}{i!}=\frac{e^x+e^{rx}+e^{r^2x}}{3}i=0∑∞[3∣i]i!xi=31i=0∑∞i!(r0+ri+r2i)xi=3ex+erx+er2x

所以答案就是(ex+erx+er2x3)k×n!(\frac{e^x+e^{rx}+e^{r^2x}}{3})^k\times n!(3ex+erx+er2x)k×n!的nnn次项系数。
将上式大力展开得：
n!×13k∑i=0k∑j=0k−iCki×Ck−ij×eix×ejrx×e(k−i−j)r2xn!\times\frac{1}{3^k}\sum_{i=0}^{k}\sum_{j=0}^{k-i}C_{k}^{i}\times C_{k-i}^{j}\times e^{ix}\times e^{jrx}\times e^{(k-i-j)r^2x}n!×3k1i=0∑kj=0∑k−iCki×Ck−ij×eix×ejrx×e(k−i−j)r2x

=n!×13k∑i=0k∑j=0k−iCki×Ck−ij×e(i+jr+(k−i−j)r2)x=n!\times\frac{1}{3^k}\sum_{i=0}^{k}\sum_{j=0}^{k-i}C_{k}^{i}\times C_{k-i}^{j}\times e^{(i+jr+(k-i-j)r^2)x}=n!×3k1i=0∑kj=0∑k−iCki×Ck−ij×e(i+jr+(k−i−j)r2)x

e(i+jr+(k−i−j)r2)xe^{(i+jr+(k-i-j)r^2)x}e(i+jr+(k−i−j)r2)x的nnn次项系数为(i+jr+(k−i−j)r2)nn!\frac{(i+jr+(k-i-j)r^2)^n}{n!}n!(i+jr+(k−i−j)r2)n

所以最终答案为：
=13k∑i=0k∑j=0k−iCki×Ck−ij×(i+jr+(k−i−j)r2)n=\frac{1}{3^k}\sum_{i=0}^{k}\sum_{j=0}^{k-i}C_{k}^{i}\times C_{k-i}^{j}\times (i+jr+(k-i-j)r^2)^n=3k1i=0∑kj=0∑k−iCki×Ck−ij×(i+jr+(k−i−j)r2)n

后话：其实，对于d=2d=2d=2的情况，–1就是二次单位根。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=19491001;
const int r=18827933;
const int K=500010;
int read(){int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;
}
int add(int a,int b){return a+b>=mod?a+b-mod:a+b;};
int dec(int a,int b){return a<b?a-b+mod:a-b;}
int mul(int a,int b){return 1ll*a*b%mod;}
int ksm(int a,int b){int res=1;while(b){if(b&1) res=mul(res,a);b>>=1;a=mul(a,a);}return res;
}
int inv(int a){return ksm(a,mod-2);
}
int n,k,d,fac[K],ifac[K];
int C(int n,int m){return mul(fac[n],mul(ifac[n-m],ifac[m]));
}
int main(){n=read();k=read();d=read();if(d==1){printf("%d\n",ksm(k,n));return 0;}fac[0]=1;for(int i=1;i<=k;i++) fac[i]=mul(fac[i-1],i);ifac[k]=inv(fac[k]);for(int i=k;i>=1;i--) ifac[i-1]=mul(ifac[i],i);if(d==2){int ans=0;for(int i=0;i<=k;i++) ans=add(ans,mul(C(k,i),ksm(2*i-k,n)));ans=mul(ans,inv(ksm(2,k)));printf("%d\n",ans);return 0;}int ans=0;for(int i=0;i<=k;i++){int sum=0;for(int j=0;j<=k-i;j++)sum=add(sum,mul(C(k-i,j),ksm(((ll)r*r*i%mod+(ll)r*j%mod+k-i-j)%mod,n)));ans=add(ans,mul(sum,C(k,i)));}ans=mul(ans,inv(ksm(3,k)));printf("%d\n",ans);return 0;
}