前置知识单位根反演自己去浅谈单位根反演看（此外可能需要一定的生成函数的姿势）
首先一看(d)这么小，那我们来分类讨论一下吧
当(d=1)时，显然答案就是(k^n)
当(d=2)时，如果你知道可重排列的指数型生成函数：
[G(x)=sum_{i=0} frac{x^{2i}}{(2i)!}
]
那么就跳过以下部分直接去看转化，我们来推导一下这个生成函数
直接搞一个DP，设(f_{i,j})表示前(i)个复读机选了(j)个时间的方案数，转移的时候枚举这个复读机复读了几次
[f_{i,j}=sum_{k=0}^j [2|k] C_{n-(j-k)}^k f_{i-1,j-k}
]
把(f_i)看做一个指数型生成函数，那么每次加入一个复读机相当于乘上了一个上面的(sum_{i=0} [2|i] frac{x^i}{i!})
由于我们这里的生成函数没有考虑每次操作间的顺序，因此乘上(n!)有：
[Ans=n!cdot G(x)^k[x^n]
]
然后这个(G(x))细看有点熟悉，稍加推导会发现它就是(frac{e^x+e^{-x}}{2})，那么这里还是直接用二项式定理展开
[Ans=n!cdot (frac{(e^x+e^{-x})^k}{2^k})=n!cdot(frac{sum_{i=0}^k C_k^i cdot e^{{2i-k}x}}{2^k})
]
然后我们求([x^n])，泰勒逆展开之后会发现分母里的(n!)和前面的(n!)抵消了，因此
[Ans=frac{sum_{i=0}^k C_k^icdot (2i-k)^n}{2^k}
]
直接(O(k))计算即可
当(d=3)时，和上面类似我们可以直接得出此时的(G(x)=sum_{i=0} [3|i] frac{x^i}{i!})，然后来一发单位根反演：
[G(x)=sum_{i=0}frac{1}{3}sum_{j=0}^2 frac{omega_3^{ij}cdot x^i}{i!}
]
[=frac{1}{3}sum_{j=0}^2 sum_{i=0}frac{(omega_3^jcdot x)^i}{i!}
]
快乐地泰勒逆展开一下就有：
[G(x)=frac{1}{3}sum_{j=0}^2 e^{omega_3^jx}
]
然后还是那个式子：
[Ans=n!cdot G(x)^k[x^n]
]
[=n!cdot (frac{sum_{a+b+c=k} (e^{(aomega_3^0+bomega_3^1+comega_3^2)x})cdot frac{k!}{a!b!c!}}{3^k})
]
还是泰勒逆展开，最后就有：
[Ans=frac{sum_{a+b+c=k} (aomega_3^0+bomega_3^1+comega_3^2)cdot frac{k!}{a!b!c!}}{3^k}
]
然后我们会发现(d=2)的情况其实也可以这么做，不过由于此时后面那个系数就是组合数，而且此时单位根就是(pm 1)，因此也可以推回上面的式子
最后说句闲话，这道题用(d=3)时的方法推一推就可以得到一个广泛的做法，此时的复杂度是(O(k^{d-1}))
打了好多公式累死了233
#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=500005,mod=19491001;
int n,k,d,fact[N],inv[N],ans,cnt,pri[100],w0,w1,w2;
inline int quick_pow(int x,int p=mod-2,int mul=1)
{
for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline void inc(int &x,CI y)
{
if ((x+=y)>=mod) x-=mod;
}
inline void init(CI n)
{
RI i; for (fact[0]=i=1;i<=n;++i) fact[i]=1LL*fact[i-1]*i%mod;
for (inv[n]=quick_pow(fact[n]),i=n-1;~i;--i) inv[i]=1LL*inv[i+1]*(i+1)%mod;
}
inline int C(CI n,CI m)
{
return 1LL*fact[n]*inv[m]%mod*inv[n-m]%mod;
}
inline int getroot(CI p)
{
RI i,j; int x=p-1; for (i=2;i*i<=x;++i) if (x%i==0)
for (pri[++cnt]=i;x%i==0;x/=i); for (i=2;;++i)
{
bool flag=1; for (j=1;j<=cnt;++j)
if (quick_pow(i,(p-1)/pri[j])==1) { flag=0; break; }
if (flag) return quick_pow(i,(mod-1)/3);
}
}
int main()
{
RI i,j; scanf("%d%d%d",&n,&k,&d); switch (d)
{
case 1:
printf("%d",quick_pow(k,n)); break;
case 2:
for (init(k),i=0;i<=k;++i) inc(ans,1LL*C(k,i)*quick_pow(2*i-k,n)%mod);
printf("%d",1LL*ans*quick_pow(2,mod-1-k)%mod); break;
case 3:
w0=1; w1=getroot(mod); w2=1LL*w1*w1%mod; init(k);
for (i=0;i<=k;++i) for (j=0;i+j<=k;++j)
inc(ans,1LL*quick_pow((1LL*w0*i%mod+1LL*w1*j%mod+1LL*w2*(k-i-j)%mod)%mod,n)*
fact[k]%mod*inv[i]%mod*inv[j]%mod*inv[k-i-j]%mod);
printf("%d",1LL*ans*quick_pow(3,mod-1-k)%mod); break;
}
return 0;
}