做法 1
可以将其看成一个 $(n+1) \times (m+1)$ 的表格，其中 $f(0,m) = 0,f(n,0) = n^K$。这个递推式的组合意义就是每次可以走向量 $(0,1)$ 和 $(1,0)$，且第一次强制走 $(0,1)$ ，走到 $(n,m)$ 的方案数，这个方案数就是 $f(n,0)$ 贡献到 $f(n,m)$ 的次数。所以可以得到答案是：
$$ \sum_{i=1}^n \binom{n-i+m-1}{m-1}i^K $$
设 $m' = m-1$，首先我们会尝试考虑将 $K$ 拆斯特林数，但是 $K$ 太大了看起来无法处理斯特林数，所以我们考虑拆前面的组合数。显然前面的组合数可以写成一个 $m'$ 次多项式 $G(n)=\binom{n+m'-1}{m'}$，现在答案也就是：
$$ \sum_{i=1}^n G(n-i)i^K \tag 1 $$
考虑继续拆式子，设 $G(x) = \sum_{i=0}^{m'} g_i x^i$，带入可以得到：
$$ \sum_{i=0}^{m'} g_i \sum_{j = 1}^n (n-j)^i j^k $$
那么我们枚举 $i$，现在重点是计算后面的东西 $F(n,k_1,k_2) = \sum_{i=1}^n (n-i)^{k_1}i^{k_2}$。由于 $k_1 \leq 30$，考虑我们对其再进行二项式展开，得到：
$$ F(n,k_1,k_2) = \sum_{i=0}^{k_1}\binom{k_1}{i} n^{i} (-1)^{k_1-i} \sum_{j=1}^n j^{k_2+k_1-i} $$
后面就是个自然数幂求和了，直接拉格朗日插值即可。复杂度 $O(mK)$。但是有个严重的问题是常数太大。
一些我使用的卡常技巧：
- 初始化拉格朗日插值的时候要求 $O(m)$ 个 $i^k$，这里可以线性筛，也提前计算好前缀和，这样可以在内存访问连续的时候计算好
- 拉格朗日插值很多次都是重复的，本质只会有 $O(m)$ 个不同的询问，可以记忆化。
- 拉格朗日插值要算个 $pre_i = \prod_{i}(x-i+1)$，由于询问来的是递增的，所以 $pre_i$ 不用每次都重新算
- 最终式子中有一项贡献是 $pre_{i-1} \times inv_{i-1}$，这个可以乘到 $pre_{i-1}$ 上，减少取模常数。
注意特判 $n=0,m=0$ 就好了。
做法 2
注意到 $(1)$ 式中 $G$ 是 $m'$ 次多项式，$i^K$ 是 $K$ 次多项式，而这个形式实际上是个卷积，所以它其实是个 $K+m'$ 次多项式，再求和大概就是 $K+m'+1$ 次多项式，直接暴力求出 $f(0,m),\ldots,f(K+m'+1,m)$ 就好了。这里其实可以直接递推！比赛的时候一直以为暴力不能求。
复杂度还是 $O(Km)$，但是常数很小。
my code:
#include <bits/stdc++.h>
#define fi first
#define se second
#define DB double
#define U unsigned
#define P std::pair
#define LL long long
#define LD long double
#define pb emplace_back
#define pf emplace_front
#define MP std::make_pair
#define SZ(x) ((int)x.size())
#define all(x) x.begin(),x.end()
#define CLR(i,a) memset(i,a,sizeof(i))
#define FOR(i,a,b) for(int i = a;i <= b;++i)
#define ROF(i,a,b) for(int i = a;i >= b;--i)
#define DEBUG(x) std::cerr << #x << '=' << x << std::endl
const int MAXN = 2500000 + 500;
const int ha = 1e9 + 7;
inline int qpow(int a,int n=ha-2){
int res = 1;
while(n){
if(n & 1) res = 1ll*res*a%ha;
a = 1ll*a*a%ha;
n >>= 1;
}
return res;
}
LL n;int m,k;
int fac[MAXN],inv[MAXN];
inline void add(int &x,int y){
x += y-ha;x += x>>31&ha;
}
int a[60];
int pw[MAXN];
int y[MAXN];
int dp[MAXN];
inline int C(int n,int m){
return n < 0 || m < 0 || n < m ? 0 : 1ll*fac[n]*inv[m]%ha*inv[n-m]%ha;
}
int pre[MAXN],suf[MAXN];
int pww[30][MAXN];
bool p[MAXN];
int prime[MAXN],cnt;
inline void prework(){
FOR(i,0,29) pww[i][1] = 1;
FOR(i,2,MAXN-1){
if(!p[i]){
prime[++cnt] = i;
pww[0][i] = qpow(i,k);
FOR(j,1,29) pww[j][i] = 1ll*pww[j-1][i]*i%ha;
}
for(int j = 1;j <= cnt && 1ll*i*prime[j] < MAXN;++j){
p[i*prime[j]] = 1;
FOR(k,0,29) pww[k][i*prime[j]] = 1ll*pww[k][i]*pww[k][prime[j]]%ha;
if(!(i%prime[j])) break;
}
}
FOR(i,2,MAXN-1) FOR(j,0,29) add(pww[j][i],pww[j][i-1]);
}
inline int gao(int m){
if(dp[m-::k]) return dp[m-::k];
int ans = 0;
FOR(i,1,m+2) y[i] = pww[m-::k][i];//qpow(i,m);
// FOR(i,2,m+2) add(y[i],y[i-1]);
int t = n%ha;
if(m == ::k){
pre[0] = suf[m+3] = 1;
FOR(i,1,m+2) pre[i] = 1ll*(t-i+ha)*pre[i-1]%ha;
ROF(i,m+2,1) suf[i] = 1ll*(t-i+ha)*suf[i+1]%ha;
FOR(i,1,m+2) pre[i] = 1ll*pre[i]*inv[i]%ha;
}
else{
pre[m+2] = 1ll*pre[m+1]*(t-(m+2)+ha)%ha*fac[m+1]%ha*inv[m+2]%ha;
suf[m+3] = 1;
ROF(i,m+2,1) suf[i] = 1ll*(t-i+ha)*suf[i+1]%ha;
}
FOR(i,1,m+2){
int c = 1ll*y[i]*pre[i-1]%ha*suf[i+1]%ha*inv[m+2-i]%ha;
if((m+2-i)&1) c = ha-c;
add(ans,c);
}
return dp[m-::k] = ans;
}
inline int calc(int k2){
int ans = 0;
FOR(i,0,k2){
int c = 1ll*C(k2,i)*qpow(n%ha,i)%ha*gao(k+k2-i)%ha;
if((k2-i)&1) c = ha-c;
add(ans,c);
}
return ans;
}
int main(){
fac[0] = 1;FOR(i,1,MAXN-1) fac[i] = 1ll*fac[i-1]*i%ha;
inv[MAXN-1] = qpow(fac[MAXN-1]);ROF(i,MAXN-2,0) inv[i] = 1ll*inv[i+1]*(i+1)%ha;
scanf("%lld%d%d",&n,&m,&k);
prework();
if(n == 0) return puts("0"),0;
if(m == 0) return printf("%d\n",qpow(n%ha,k)),0;
a[0] = 1;
FOR(i,0,m-2){// * (x-i)
ROF(j,59,1) a[j] = (a[j-1]-1ll*i*a[j]%ha+ha)%ha;
a[0] = (ha-1ll*a[0]*i%ha)%ha;
}
FOR(i,0,m-1) a[i] = 1ll*a[i]*inv[m-1]%ha;
int ans = 0;
// FOR(i,1,n) add(ans,1ll*qpow(i,k)*C(n-i+m-1,m-1)%ha);
FOR(i,0,m-1){
int sm = 0;
FOR(j,i,m-1) add(sm,1ll*a[j]*qpow(m-1,j-i)%ha*C(j,i)%ha);
a[i] = sm;
}
FOR(i,0,m-1){
add(ans,1ll*a[i]*calc(i)%ha);
}
printf("%d\n",ans);
return 0;
}
