直接摘抄 OI-wiki 了。
第二类斯特林数
第二类斯特林数（斯特林子集数）\(\begin{Bmatrix}n\\ k\end{Bmatrix}\)，也可记做 \(S(n,k)\)，表示将 \(n\) 个两两不同的元素，划分为 \(k\) 个互不区分的非空子集的方案数。
递推式
\[\begin{Bmatrix}n\\ k\end{Bmatrix}=\begin{Bmatrix}n-1\\ k-1\end{Bmatrix}+k\begin{Bmatrix}n-1\\ k\end{Bmatrix} \]边界是 \(\begin{Bmatrix}n\\ 0\end{Bmatrix}=[n=0]\)。
考虑用组合意义来证明。
我们插入一个新元素时，有两种方案：
- 将新元素单独放入一个子集，有 \(\begin{Bmatrix}n-1\\ k-1\end{Bmatrix}\) 种方案；
- 将新元素放入一个现有的非空子集，有 \(k\begin{Bmatrix}n-1\\ k\end{Bmatrix}\) 种方案。
根据加法原理，将两式相加即可得到递推式。
普通幂转下降幂
我们记下降阶乘幂 \(x^{\underline{n}}=\dfrac{x!}{(x-n)!}=\prod_{k=0}^{n-1} (x-k)\)。
则可以利用下面的恒等式将普通幂转化为下降幂：
\[x^n=\sum_{k} \begin{Bmatrix}n\\ k\end{Bmatrix} x^{\underline{k}} \]含义是：有 \(n\) 个有标号小球，想要放进 \(x\) 个有标号盒子中，问有多少种方案。第一种算法是 \(x^n\)。第二种算法则考虑最终有 \(k\) 个盒子是非空的，小球划分为 \(k\) 个无标号集合，从盒子中有标号地选 \(k\) 个盒子出来。
吸收恒等式
\[m^{\underline k}\binom n m=n^{\underline k}\binom {n-k}{m-k} \]有 \(n\) 个有标号小球，先无序地选出 \(m\) 个，再从 \(m\) 个中有序地选出 \(k\) 个，等价于选从 \(n\) 个中有序地选出 \(k\) 个，再从剩下的小球中无序地选出 \(m-k\) 个。
题目描述
众所周知，小葱同学擅长计算，尤其擅长计算组合数。小葱现在希望你计算
\[\left(\sum_{k=0}^{n}f(k)\times x^k\times \binom{n}{k}\right)\bmod p \]的值。其中 \(n\), \(x\), \(p\) 为给定的整数，\(f(k)\) 为给定的一个 \(m\) 次多项式 \(f(k) = a_0 + a_1k + a_2k^2 + \cdots + a_mk^m\)。\(\binom{n}{k}\) 为组合数，其值为 \(\binom{n}{k} = \frac{n!}{k!(n-k)!}\)。
对于所有测试数据：\(1\le n, x, p \le 10^9, 0\le a_i\le 10^9, 0\le m \le \min(n,1000)\)。
solution
\(f(k)\) 转为下降幂形式 \(\sum_{i=0}^m b_ik^{\underline i}\)。则有
\[\begin{aligned} &=\sum_{k=0}^n\sum_{i=0}^mb_ik^{\underline i}x^k\binom n k\\ &=\sum_{k=0}^n\sum_{i=0}^mb_in^{\underline i}x^k\binom {n-i} {k-i}\\ &=\sum_{i=0}^mb_in^{\underline i}\sum_{k=0}^nx^k\binom {n-i} {k-i}\\ &=\sum_{i=0}^mb_in^{\underline i}x^i\sum_{k=0}^{n-i}x^k\binom {n-i} {k}\\ &=\sum_{i=0}^mb_in^{\underline i}x^i(1+x)^{n-i} \end{aligned} \]足以计算。依实现，可以 \(O(m^2\log n)\) 或 \(O(m^2)\)。
code
#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
using LL = long long;
template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <int id>
struct modint {/*{{{*/
static int mod;
static unsigned umod;
static void setmod(int p) { mod = umod = p; }
unsigned v;
modint() : v(0) {}
template <class T, must_int<T> = 0>
modint(T x) {
x %= mod;
v = x < 0 ? x + mod : x;
}
modint operator+() const { return *this; }
modint operator-() const { return modint() - *this; }
friend int raw(const modint &self) { return self.v; }
friend ostream& operator<<(ostream& os, const modint &self) {
return os << raw(self);
}
modint &operator+=(const modint &rhs) {
v += rhs.v;
if (v >= umod) v -= umod;
return *this;
}
modint &operator-=(const modint &rhs) {
v -= rhs.v;
if (v >= umod) v += umod;
return *this;
}
modint &operator*=(const modint &rhs) {
v = 1ull * v * rhs.v % umod;
return *this;
}
modint &operator/=(const modint &rhs) {
assert(rhs.v);
return *this *= qpow(rhs, mod - 2);
}
template <class T, must_int<T> = 0>
friend modint qpow(modint a, T b) {
modint r = 1;
for (; b; b >>= 1, a *= a)
if (b & 1) r *= a;
return r;
}
friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
bool operator==(const modint &rhs) const { return v == rhs.v; }
bool operator!=(const modint &rhs) const { return v != rhs.v; }
};/*}}}*/
template <int id>
unsigned modint<id>::umod;
template <int id>
int modint<id>::mod;
using mint = modint<-1>;
int n, m;
mint X, b[1010], S[1010][1010];
int main() {
#ifndef LOCAL
cin.tie(nullptr)->sync_with_stdio(false);
#endif
cin >> n >> X.v >> mint::mod >> m;
mint::setmod(mint::mod);
X.v %= mint::mod;
S[0][0] = 1;
for (int i = 1; i <= m; i++) {
S[i][0] = 0;
for (int j = 1; j <= i; j++) S[i][j] = S[i - 1][j - 1] + S[i - 1][j] * j;
}
for (int i = 0; i <= m; i++) {
int x;
cin >> x;
for (int j = 0; j <= i; j++) b[j] += S[i][j] * x;
}
mint ans = 0;
mint now = 1;
for (int i = 0; i <= m; i++) {
ans += b[i] * now * qpow(X, i) * qpow(1 + X, n - i);
now *= n - i;
}
cout << ans << endl;
return 0;
}