Link
考虑单个元素的EGF，显然为(f(x)=sumlimits_{ige0}[d|i]frac{x^i}{i!}=frac1dsumlimits_{i=0}^{d-1}exp(omega_d^ix))。
那么答案就是([x^n]f(x)^k)。
考虑分圆多项式(Phi_d(x)=prodlimits_{kin[0,d)wedgegcd(d,k)=1}(x-omega_d^k)=prodlimits_{k|d}(x^{frac dk}-1)^{mu(k)})。
那么(omega_d^k=(x^kmodPhi_d(x))mid_{x=omega})。
很显然(deg(Phi_d(x))=varphi(d))，且(Phi_d(x))的系数都是整数，因此(omega_d^k)可以由(omega_d^0,cdots,omega_d^{varphi(d)-1})以整数线性表示。
实际上可以证明，对于域(mathbb N)上的线性空间(V=operatorname{span}(omega_d^0,cdots,omega_d^{d-1}))一定有(dim V=varphi(d))。
这样我们就可以把(f(x))写成(F(exp x,exp(omega x),cdots,exp(omega_d^{d-1}x)))的形式了。
设(G=F^k)，那么(Ffrac{partial G}{partial x}=kGfrac{partial F}{partial x})。
由这个式子列出递推式然后解得系数是trivial的。
时间复杂度为(O(k^{varphi(d)}log n))，(log n)是快速幂的。
当(d=2)时可以通过线性筛幂函数做到(O(k))。
#include<cstdio>
using i64=long long;
const int N=2007,P=19491001;
int n,k,d;i64 w,ans,fac[N],ifac[N],inv[N];
int abs(int x){return x<0? -x:x;}
void inc(i64&a,i64 b){a+=b-P,a+=a>>63&P;}
i64 pow(i64 a,i64 b){i64 c=1;for(a%=P;b;b>>=1,a=a*a%P)if(b&1)c=c*a%P;return c;}
i64 binom(int n,int m){return fac[n]*ifac[m]%P*ifac[n-m]%P;}
void init(int n=2000)
{
fac[0]=ifac[0]=inv[0]=fac[1]=ifac[1]=inv[1]=1,w=pow(7,(P-1)/d);
for(int i=2;i<=n;++i) fac[i]=i*fac[i-1]%P,inv[i]=(P-P/i)*inv[P%i]%P,ifac[i]=inv[i]*ifac[i-1]%P;
}
namespace Sub1
{
void solve()
{
for(int a=0;a<=k;++a) for(int b=2-((k^a)&1);a+b<=k;b+=2) inc(ans,pow(a*w+b,n)*binom(k,(k+a+b)/2)%P*binom(k,(k+abs(a-b))/2)%P);
printf("%lld",4*pow(d,P-k-1)*ans%P);
}
}
namespace Sub2
{
i64 a[N][N],b[N][N];
void solve()
{
for(int i=0;i<=k;++i) a[0][i]=a[i][0]=binom(k,i),b[i][0]=i*a[i][0]%P;
for(int i=1;i<=k+k;++i)
for(int j=1;j<=i&&j<=k;++j)
{
i64 x=k*a[i-1][j]-b[i-1][j]-b[i][j-1];
if(j>=2) x+=k*a[i-1][j-2]-b[i-1][j-2];
if(i>=2) x+=2*k*a[i-2][j-1]-b[i-2][j-1];
if(i>=2&&j>=2) x+=2*k*a[i-2][j-2]-b[i-2][j-2];
x=(x+4*P)%P,b[i][j]=x,a[i][j]=inv[i]*x%P;
if(i<=k) a[j][i]=a[i][j],b[j][i]=j*a[j][i]%P;
}
for(int i=0;i<=k+k;++i) for(int j=0;j<=k;++j) inc(ans,(j==k? 1:2)*pow((P+i-k)*w+k-j,n)*a[i][j]%P);
printf("%lld",ans*pow(d,P-k-1)%P);
}
}
int main()
{
scanf("%d%d%d",&n,&k,&d),init();
if(n%d) puts("0"); else if(d==4) Sub1::solve(); else if(d==6) Sub2::solve();
}