从阶乘幂到斯特林数 - caijianhong - 博客园 (cnblogs.com)
题目描述
Crash 小朋友最近迷上了一款游戏——文明5 (Civilization V)。在这个游戏中，玩家可以建立和发展自己的国家，通过外交和别的国家交流，或是通过战争征服别的国家。
现在 Crash 已经拥有了一个 \(n\) 个城市的国家，这些城市之间通过道路相连。由于建设道路是有花费的，因此 Crash 只修建了 \(n-1\) 条道路连接这些城市，不过可以保证任意两个城市都有路径相通。
在游戏中，Crash 需要选择一个城市作为他的国家的首都，选择首都需要考虑很多指标，有一个指标是这样的：
\[S(i) = \sum_{j = 1}^{n}{\rm dist}(i, j) ^ k \]其中 \(S(i)\) 表示第 \(i\) 个城市的指标值，\({\rm dist}(i, j)\) 表示第 \(i\) 个城市到第 \(j\) 个城市需要经过的道路条数的最小值，\(k\) 为一个常数且为正整数。
因此 Crash 交给你一个简单的任务：给出城市之间的道路，对于每个城市，输出这个城市的指标值，由于指标值可能会很大，所以你只需要输出这个数 \(\bmod\ 10007\) 的值。
对于 \(100 \%\) 的数据，\(1\le n\le 5\times 10^4\)，\(1\le k\le 150\)。
solution
为什么这东西不能直接换根 dp？你考虑你想要的操作无非就是：给你一个集合 \(S\)，用
\[\sum_{x\in S}x^k \]求出
\[\sum_{x\in S}(x+1)^k \]然而假如你只管 \(k\) 次方的信息，你减一下就会发现
\[(x+1)^k-x^k=\sum_{i=0}^{k-1}\binom k i x^i \]这个东西是二项式定理形式，展开以后有 \(k\) 项，都比 \(k\) 次要低。这 \(k\) 项你都要维护，维护一项需要 \(O(k)\)，合着你需要总共 \(O(k^2)\) 的时间维护信息，过不去。
\((x+1)^k-x^k\) 是一个形如差分的东西，回想起《具体数学》第二章中介绍的有限微积分，其定义
\[\Delta f(x)=f(x+1)-f(x) \]并发现 \(\Delta x^k\) 没有很好的形式，但是 \(\Delta x^\underline k\) 有
\[\Delta x^\underline k=(x+1-(x-k+1))x^\underline {k-1}=kx^\underline {k-1} \]因为我们熟知下降幂与普通幂的转换
\[x^n=\sum_{k}\begin{Bmatrix}n\\ k\end{Bmatrix}x^\underline k \]加一个求和符号不会影响什么东西
\[\sum_{x\in S}x^n=\sum_{k}\begin{Bmatrix}n\\ k\end{Bmatrix}\sum_{x\in S}x^\underline k \]那我们就维护 \(\sum_{x\in S}x^\underline k\) 好了，就是所有 \(\leq k\) 次下降阶乘幂都维护一下，这样就能实现我们想要的东西了
\[\sum_{x\in S}(x+1)^\underline k=k\sum_{x\in S}x^\underline {k-1}+\sum_{x\in S}x^\underline k \]可以 \(O(k)\) 维护出集合中所有数字加一后的信息，可以换根 dp。
\(f_{u, k}\) 维护了 \(\sum_{x\in subtree(u)}x^\underline k\)。记 \(A(f_u)\) 表示将所维护的所有 \(x\) 加一。
转移：\(f_u=\sum_v A(f_v)\)，再加上自己的一个 \(0\)（意思是 \(0^\underline k=[k=0]\)）。
换根：假设 \(f_u\) 现在拿着是全局的信息，那么将 \(A(f_u-A(f_v))\) 加到 \(f_v\) 上，然后递归 \(v\)。
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 <unsigned umod>
struct modint {/*{{{*/
static constexpr int mod = umod;
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 = 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; }
};/*}}}*/
using mint = modint<10007>;
int n, k;
mint S[210][210];
valarray<mint> a[50010];
vector<int> g[50010];
valarray<mint> add1(valarray<mint> f) {
for (int i = k; i >= 1; i--) f[i] += f[i - 1] * i;
return f;
}
mint getValue(const valarray<mint>& f) {
mint ret = 0;
for (int i = 0; i <= k; i++) ret += S[k][i] * f[i];
return ret;
}
void dfs(int u, int fa) {
a[u][0] = 1;
for (int v : g[u]) if (v != fa) {
dfs(v, u);
a[u] += add1(a[v]);
}
}
void chr(int u, int fa) {
for (int v : g[u]) if (v != fa) {
a[v] += add1(a[u] - add1(a[v]));
chr(v, u);
}
}
int main() {
#ifndef LOCAL
cin.tie(nullptr)->sync_with_stdio(false);
#endif
cin >> n >> k;
S[0][0] = 1;
for (int i = 1; i <= k; i++) {
for (int j = 1; j <= i; j++) S[i][j] = S[i - 1][j] * j + S[i - 1][j - 1];
}
for (int i = 1; i <= n; i++) a[i].resize(k + 1);
for (int i = 1, u, v; i < n; i++) cin >> u >> v, g[u].push_back(v), g[v].push_back(u);
dfs(1, 0);
chr(1, 0);
for (int i = 1; i <= n; i++) cout << getValue(a[i]) << endl;
return 0;
}