拉格朗日插值
拉格朗日插值
一般性插值
给出 \(n\) 个不同的点值 \((x_{0 \sim n} , y_{0 \sim n})\) 和常数 \(k\) ，求一个 \(n - 1\) 次多项式 \(f\) ，满足其经过给出的 \(n\) 个点，并输出 \(f(k)\) 。
\(n \leq 2000\)
对于第 \(i\) 个点，考虑构造一个多项式 \(F_i(k) = \begin{cases} y_i & k = x_i \\ 0 & k \not = x_i \end{cases}\) ，则 \(F(k) = \sum_{i = 0}^n F_i(k)\) 即为所求。
可以构造：
\[F_i(k) = y_i \prod_{j \ne i} \frac{k - x_j}{x_i - x_j}
\]
解释：当 \(k = x_j\) 时，分子的累乘肯定有一项为 \(0\) ，反之分子分母累乘后肯定都会被约成 \(1\) 。
于是得到：
\[F(k) = \sum_{i = 0}^n y_i \prod_{j \ne i} \frac{k - x_j}{x_i - x_j}
\]
忽略求逆元的复杂度，时间复杂度 \(O(n^2)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353;
const int N = 2e3 + 7;
int x[N], y[N];
int n, k;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
signed main() {
scanf("%d%d", &n, &k);
for (int i = 1; i <= n; ++i)
scanf("%d%d", x + i, y + i);
int ans = 0;
for (int i = 1; i <= n; ++i) {
int res = y[i];
for (int j = 1; j <= n; ++j)
if (j != i)
res = 1ll * res * dec(k, x[j]) % Mod * mi(dec(x[i], x[j]), Mod - 2) % Mod;
ans = add(ans, res);
}
printf("%d", ans);
return 0;
}
横坐标连续的情况
把 \(x_i\) 替换为 \(i\) ，新的式子即为：
\[F(k) = \sum_{i = 0}^n y_i \prod_{i \neq j} \frac{k - j}{i - j}
\]
考虑优化。观察到分子等于 \(\frac{\prod_{j = 0}^n (k - j)}{k - i}\) ，分母等于 \((-1)^{n - i} (i - 1)! (n - i)!\) ，于是得到：
\[G(k) = \sum_{i = 0}^n y_i \dfrac{\prod_{j = 0}^n (k - j)}{(k - i) (-1)^{n - i} (i - 1)! (n - i)!}
\]
重心拉格朗日插值
观察式子：
\[F(k) = \sum_{i = 0}^n y_i \prod_{i \not = j} \dfrac{k - x_j}{x_i - x_j}
\]
令 \(g = \prod_{i = 0}^n (k - x_i)\) ，则：
\[G(k) = g \sum_{i = 0}^n \prod_{i \not = j} \dfrac{y_i}{(k - x_i)(x_i - x_j)}
\]
再令 \(t_i = \frac{y_i}{\prod_{j \not = i} x_i - x_j}\) ，则：
\[G(k) = g \sum_{i = 0}^n \dfrac{t_i}{k - x_i}
\]
这样每次新加入一个点的时候只需要 \(O(n)\) 计算 \(g, t_i\) 即可。
反推系数
某一类题目给出了多项式的若干点值，需要求出每一项的系数。
先观察拉插的式子：
\[\begin{aligned}
F(x) &= \sum_{i = 0}^n y_i \prod_{i \not = j} \frac{x - x_j}{x_i - x_j} \\
&= \sum_{i = 0}^n \frac{y_i}{\prod_{i \not = j} (x_i - x_j)} \times \frac{\prod_j (x - x_j)}{x - x_i}
\end{aligned}
\]
前面一项是可以暴力求的。对于后面一项，考虑先 \(O(n^2)\) 用背包 DP 求出 \(G(x) = \prod_{i = 1}^n (x - x_i)\) 的每一项系数，则后者等于 \(\frac{G(x)}{x - x_i}\) ，直接手动模拟长除法即可。总时间复杂度 \(O(n^2)\) 。
inline vector<int> Lagrange(vector<int> X, vector<int> Y, int n) {
vector<int> g(n + 1);
g[0] = 1;
for (int i = 0; i < n; ++i)
for (int j = i + 1; ~j; --j)
g[j] = add(1ll * dec(0, X[i]) * g[j] % Mod, j ? g[j - 1] : 0);
vector<int> f(n);
for (int i = 0; i < n; ++i) {
int mul = 1;
for (int j = 0; j < n; ++j)
if (i != j)
mul = 1ll * mul * dec(X[i], X[j]) % Mod;
mul = 1ll * Y[i] * mi(mul, Mod - 2) % Mod;
for (int j = n - 1, res = g[n]; ~j; --j)
f[j] = add(f[j], 1ll * mul * res % Mod), res = add(g[j], 1ll * X[i] * res % Mod);
}
return f;
}
应用
CF622F The Sum of the k-th Powers
给定 \(n, k\) ，求 \((\sum_{i = 1}^n i^k) \bmod (10^9 + 7)\) 。
\(n \leq 10^9\) ，\(k \leq 10^6\)
结论：\(\sum_{i = 1}^n i^k\) 是一个关于 \(n\) 的 \(k + 1\) 次多项式。
证明：定义数列 \(\{ a_n \}\) 的 \(p\) 阶差分数列为 \(\Delta^p f(x) = \Delta(\Delta^{p - 1} f(x))\) 。若数列 \(\{ a_n \}\) 的 \(p\) 阶阶差数列是一个非 \(0\) 的常数列，则称其为 \(p\) 阶等差数列。
定理：数列 \(\{a_n\}\) 是一个 \(p\) 阶等差数列等价于其通项公式是一个关于 \(n\) 的 \(p\) 次多项式。
证明：考虑证明其充分性，必要性的证法也是类似的。
设数列的通项公式 \(f(x) = \sum_{i = 0}^n a_i x^i\) ，则其一阶差分数列的通项公式为：
\[\begin{aligned} \Delta f(x) &= f(x + 1) - f(x) \\ &= \sum_{i = 0}^n a_i (x + 1)^i - \sum_{i = 0}^n a_i x^i \\ &= \sum_{i = 0}^n a_i [(x + 1)^i - x^i] \end{aligned} \]只考虑 \(x_n\) 这一项：
\[a_n[(x + 1)^n - x^n] \]把 \((x + 1)^n\) 用二项式定理展开，因为只考虑 \(n\) 次项，所以我们只保留 \(x^n\) ，发现其系数为 \(1\) ，刚好和后面的项消掉。
因此将数列差分一次后就会消掉最高次项，于是 \(p\) 次差分后就会消掉 \(p\) 次项，因此得证。
设数列 \(\{ a_n \}\) ，其通项公式为 \(a_n = \sum_{i = 1}^n i^k\) ，其差分数列的通项公式为 \(f(x) = x^k\) 是一个 \(k\) 次多项式，所以这个数列是一个关于 \(n\) 的 \(k + 1\) 次多项式。
带入 \(1 \sim k + 2\) 的点值，利用线性筛预处理 \(i^k\) 即可做到 \(O(n)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1e6 + 7;
int pri[N], f[N], inv[N], invfac[N], pre[N], suf[N];
bool isp[N];
int n, k, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline int sgn(int n) {
return n & 1 ? Mod - 1 : 1;
}
inline void prework(int n) {
memset(isp + 1, true, sizeof(bool) * n);
isp[1] = false, f[1] = 1;
for (int i = 2; i <= n; ++i) {
if (isp[i])
pri[++pcnt] = i, f[i] = mi(i, k);
for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
isp[i * pri[j]] = false, f[i * pri[j]] = 1ll * f[i] * f[pri[j]] % Mod;
if (!(i % pri[j]))
break;
}
}
for (int i = 2; i <= n; ++i)
f[i] = add(f[i], f[i - 1]);
inv[0] = inv[1] = 1;
invfac[0] = invfac[1] = 1;
for (int i = 2; i <= n; ++i) {
inv[i] = 1ll * (Mod - Mod / i) * inv[Mod % i] % Mod;
invfac[i] = 1ll * invfac[i - 1] * inv[i] % Mod;
}
}
signed main() {
scanf("%d%d", &n, &k);
prework(k + 2);
if (n <= k + 2)
return printf("%d", f[n]), 0;
pre[0] = 1;
for (int i = 1; i <= k + 2; ++i)
pre[i] = 1ll * pre[i - 1] * (n - i) % Mod;
suf[k + 3] = 1;
for (int i = k + 2; i; --i)
suf[i] = 1ll * suf[i + 1] * (n - i) % Mod;
int ans = 0;
for (int i = 1; i <= k + 2; ++i)
ans = add(ans, 1ll * f[i] * pre[i - 1] % Mod * suf[i + 1] % Mod *
invfac[i - 1] % Mod * sgn(k + 2 - i) % Mod * invfac[k + 2 - i] % Mod);
printf("%d", ans);
return 0;
}
P5667 拉格朗日插值2
给定一个不超过 \(n\) 次的多项式的 \(n + 1\) 个点值 \(f(0), f(1), \cdots, f(n)\)，和一个正整数 \(m\)，求 \(f(m), f(m + 1), \cdots, f(m + n)\) 。答案对 \(998244353\) 取模。
\(n \leq 1.6 \times 10^5\) ，\(n < m \leq 10^8\)
由插值公式可得：
\[\begin{aligned}
f(m + x) &= \sum_{i = 0}^n f(i) \prod_{j \not = i} \dfrac{m + x - j}{i - j} \\
&= \sum_{i = 0}^n \frac{f(i) \times \frac{(m + x)!}{(m + x - n - 1)!}}{(m + x - i) (-1)^{n - i} i! (n - i)!} \\
&= \frac{(m + x)!}{(m + x - n - 1)!} \sum_{i = 0}^n \frac{f(i)}{(m + x - i) (-1)^{n - i} i! (n - i)!}
\end{aligned}
\]
令：
\[\begin{aligned}
u_i &=
\begin{cases}
\frac{f(i)}{(-1)^{n - i} i! (n - i)!} & i \leq n \\
0 & i > n
\end{cases} \\
v_i &= \frac{1}{m - n + i}
\end{aligned}
\]
可得：
\[f(m + x) = \frac{(m + x)!}{(m + x - n - 1)!} \times (u \times v)[n + x]
\]
NTT 即可做到 \(O(n \log n)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353, rt = 3, invrt = 332748118;
const int N = 2e6 + 7;
int fac[N], inv[N], invfac[N];
int f[N], u[N], v[N];
int n, m;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline int sgn(int n) {
return n & 1 ? Mod - 1 : 1;
}
inline void prework(int n) {
fac[0] = fac[1] = 1;
inv[0] = inv[1] = 1;
invfac[0] = invfac[1] = 1;
for (int i = 2; i <= n; ++i) {
fac[i] = 1ll * fac[i - 1] * i % Mod;
inv[i] = 1ll * (Mod - Mod / i) * inv[Mod % i] % Mod;
invfac[i] = 1ll * invfac[i - 1] * inv[i] % Mod;
}
}
namespace Poly {
#define cpy(f, g, n) memcpy(f, g, sizeof(int) * n)
#define clr(f, n) memset(f, 0, sizeof(int) * n)
int rev[N];
inline int calc(int n) {
int len = 1;
while (len < n)
len <<= 1;
for (int i = 0; i < len; ++i)
rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? len >> 1 : 0);
return len;
}
inline void NTT(int *f, int n, int op) {
for (int i = 0; i < n; ++i)
if (i < rev[i])
swap(f[i], f[rev[i]]);
for (int k = 1; k < n; k <<= 1) {
int tG = mi(op == 1 ? rt : invrt, (Mod - 1) / (k << 1));
for (int i = 0; i < n; i += k << 1) {
int buf = 1;
for (int j = 0; j < k; ++j) {
int fl = f[i + j], fr = 1ll * buf * f[i + j + k] % Mod;
f[i + j] = add(fl, fr), f[i + j + k] = dec(fl, fr);
buf = 1ll * buf * tG % Mod;
}
}
}
if (op == -1) {
int invn = mi(n, Mod - 2);
for (int i = 0; i < n; ++i)
f[i] = 1ll * invn * f[i] % Mod;
}
}
inline void Mul(int *f, int n, int *g, int m, int *res) {
static int a[N], b[N];
int len = calc(n + m - 1);
cpy(a, f, n), clr(f + n, len - n);
cpy(b, g, m), clr(b + m, len - m);
NTT(a, len, 1), NTT(b, len, 1);
for (int i = 0; i < len; ++i)
a[i] = 1ll * a[i] * b[i] % Mod;
NTT(a, len, -1), cpy(res, a, n + m - 1);
}
#undef cpy
#undef clr
} // namespace Poly
signed main() {
scanf("%d%d", &n, &m);
for (int i = 0; i <= n; ++i)
scanf("%d", f + i);
prework(n * 2);
for (int i = 0; i <= n * 2; ++i) {
u[i] = (i <= n ? 1ll * f[i] * sgn(n - i) % Mod * invfac[i] % Mod * invfac[n - i] % Mod : 0);
v[i] = mi(m - n + i, Mod - 2);
}
Poly::Mul(u, n * 2 + 1, v, n * 2 + 1, u);
int mul = 1;
for (int i = m - n; i <= m; ++i)
mul = 1ll * mul * i % Mod;
for (int i = 0; i <= n; ++i) {
printf("%lld ", 1ll * mul * u[n + i] % Mod);
mul = 1ll * mul * mi(m + i - n, Mod - 2) % Mod * (m + i + 1) % Mod;
}
return 0;
}
[AGC063E] Child to Parent
给定一棵树，根为 \(1\) ，第 \(i\) 个点的权值为 \(a_i\) 。一次操作可以选择一个 \(u\) 满足 \(u \neq 1\) 且 \(a_u > 0\) ，然后令 \(a_i \to a_i - 1\) ，\(a_{fa_i} \to a_{fa_i} + k\) 。
求若干次操作后不同局面的数量，两个局面不同当且仅当存在一个 \(a_i\) 不同。
\(n \leq 300\)
设 \(c_i\) 表示 \(i\) 位置的操作次数，由于 \(a\) 和 \(c\) 构成双射，因此只需考虑对 \(c\) 计数。对于每个 \(u\) ，限制为 \(0 \le c_i \le a_u + k \times \sum_{v \in son(u)} c_v\) 。
设 \(f_{u, i}\) 表示 \(u\) 子树所有方案是上传权值的 \(i\) 次幂之和。先考虑 \(f_{u, 0}\) 的求法，拆开每个子树的贡献得到：
\[f_{u, 0} = (1 + a_u) \times \left( \prod_{v \in son(u)} f_{v, 0} \right) + \sum_{v \in son(u)} k \times f_{v, 1} \times \prod_{w \in son(u), w \neq v} f_{w, 0}
\]
设 \(c_i\) 表示 \(i\) 位置的操作次数，由于 \(a\) 和 \(c\) 构成双射，因此只需考虑对 \(c\) 计数。对于每个 \(u\) ，限制为 \(0 \le c_i \le a_u + k \times \sum_{v \in son(u)} c_v\) 。
设 \(f_{u, i}\) 表示 \(u\) 子树所有方案是上传权值的 \(i\) 次幂之和。先考虑 \(f_{u, 0}\) 的求法，拆开每个子树的贡献得到：
\[f_{u, 0} = (1 + a_u) \times \left( \prod_{v \in son(u)} f_{v, 0} \right) + \sum_{v \in son(u)} f_{v, 1} \times k \times \prod_{w \in son(u), w \neq v} f_{w, 0}
\]
发现要知道 \(f_{u, 0}\) 就需要知道 \(f_{v, 1}\) ，以此类推。
记 \(F_n(x) = \sum_{i = 0}^x i^n\) ，规定 \(0^0 = 1\) ，则：
\[f_{u, i} = \sum_{所有方案} F_i(1 + a_u + k \times \small\sum_{v \in son(u)} c_v)
\]
由于 \(F_n(x)\) 是一个关于 \(x\) 的 \(n + 1\) 次多项式，因此考虑维护：
\[g_{u, i} = \sum_{所有方案} (1 + a_u + k \times \small\sum_{v \in son(u)} c_v)^i
\]
即多项式中每一个 \(i\) 次幂项的和，直接拉格朗日插值即可得到 \(f\) 各项系数，那么求出 \(g\) 之后直接通过系数算出 \(f\) 即可。
求 \(g\) 的系数需要合并子树信息，直接用二项式定理即可。
时间复杂度 \(O(n^3)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353;
const int N = 3e2 + 7;
struct Graph {
vector<int> e[N];
inline void insert(int u, int v) {
e[u].emplace_back(v);
}
} G;
int fa[N], a[N], C[N][N], pw[N], f[N][N], g[N][N], h[N];
int n, m;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void prework() {
for (int i = C[0][0] = 1; i <= n; ++i)
for (int j = C[i][0] = 1; j <= i; ++j)
C[i][j] = add(C[i - 1][j], C[i - 1][j - 1]);
pw[0] = 1;
for (int i = 1; i <= n; ++i)
pw[i] = 1ll * pw[i - 1] * m % Mod;
}
namespace Lagrange {
int h[N][N], x[N], y[N], g[N];
inline void solve(int *x, int *y, int n, int *f) {
memset(g, 0, sizeof(int) * (n + 1)), g[0] = 1;
for (int i = 1; i <= n; ++i)
for (int j = i; ~j; --j)
g[j] = add(1ll * dec(0, x[i]) * g[j] % Mod, j ? g[j - 1] : 0);
memset(f, 0, sizeof(int) * n);
for (int i = 1; i <= n; ++i) {
int mul = 1;
for (int j = 1; j <= n; ++j)
if (j != i)
mul = 1ll * mul * dec(x[i], x[j]) % Mod;
mul = 1ll * y[i] * mi(mul, Mod - 2) % Mod;
for (int j = n - 1, res = g[n]; ~j; --j)
f[j] = add(f[j], 1ll * mul * res % Mod), res = add(g[j], 1ll * x[i] * res % Mod);
}
}
inline void prework() {
for (int i = 0; i <= n; ++i) {
y[0] = !i;
for (int j = 1; j <= i + 2; ++j)
x[j] = j, y[j] = add(y[j - 1], mi(j, i));
solve(x, y, i + 2, h[i]);
}
}
inline int query(int x, int n) {
int ans = 0;
for (int i = 0, pw = 1; i <= n + 1; ++i, pw = 1ll * pw * x % Mod)
ans = add(ans, 1ll * h[n][i] * pw % Mod);
return ans;
}
} // namespace Lagrange
void dfs(int u) {
if (G.e[u].empty()) {
for (int i = 0; i <= n; ++i)
f[u][i] = Lagrange::query(a[u], i);
return;
}
for (int v : G.e[u])
dfs(v);
if (u == 1)
return;
for (int i = 0; i <= n; ++i)
g[u][i] = mi(a[u], i);
for (int v : G.e[u]) {
memset(h, 0, sizeof(int) * (n + 1));
for (int i = 0; i <= n; ++i)
for (int j = 0; i + j <= n; ++j)
h[i + j] = add(h[i + j], 1ll * C[i + j][i] * pw[i] % Mod * f[v][i] % Mod * g[u][j] % Mod);
memcpy(g[u], h, sizeof(int) * (n + 1));
}
for (int i = 0; i <= n; ++i)
for (int j = 0; j <= i + 1; ++j)
f[u][i] = add(f[u][i], 1ll * Lagrange::h[i][j] * g[u][j] % Mod);
}
signed main() {
scanf("%d", &n);
for (int i = 2; i <= n; ++i)
scanf("%d", fa + i), G.insert(fa[i], i);
scanf("%d", &m);
for (int i = 1; i <= n; ++i)
scanf("%d", a + i);
prework(), Lagrange::prework(), dfs(1);
int ans = 1;
for (int x : G.e[1])
ans = 1ll * ans * f[x][0] % Mod;
printf("%d", ans);
return 0;
}
CF1874E Jellyfish and Hack
给出函数 \(fun\) 的伪代码：
function fun(A) if A.length > 0 let L[1 ... L.length] and R[1 ... R.length] be new arrays L.length = R.length = 0 for i = 2 to A.length if A[i] < A[1] L.length = L.length + 1 L[L.length] = A[i] else R.length = R.length + 1 R[R.length] = A[i] return A.length + fun(L) + fun(R) else return 0给定 \(n, lim\) ，求有多少 \(1 \sim n\) 的排列 \(p\) 满足 \(fun(p) \geq lim\) ，答案对 \(10^9 + 7\) 取模。
\(n \leq 200\) ，\(lim \leq 10^9\)
考虑 DP，设 \(f_{i, j}\) 表示长度为 \(i\) 的排列答案为 \(j\) 的方案数，则：
\[f_{i, j} = \sum_{k = 1}^i \binom{i - 1}{k - 1} \sum_{l = 0}^{j - i} f_{k - 1, l} \times f_{i - k, j - i - l}
\]
由于第二维是 \(O(n^2)\) 的，直接做可以做到 \(O(n^6)\) ，无法通过。
发现第二维是一个卷积的形式，考虑 OGF，令：
\[F_i(x) = \sum_{j = 0}^{\frac{n(n + 1)}{2}} f_{i, j} x^j
\]
则：
\[F_i = \sum_{j = 1}^i \binom{i - 1}{j - 1} x^i F_{j - 1} F_{i - j}
\]
由于 \(\deg F_i \leq \frac{n(n + 1)}{2}\) ，因此只要知道 \(O(n^2)\) 个位置的点值，即可利用拉格朗日插值反推出 \(F_i\) 的各项系数，每次直接平方推系数即可做到 \(O(n^4)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 2e2 + 7, M = N * N / 2;
int C[N][N];
int f[M], g[M], y[M];
int n, lim, m;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline int calc(int x) {
static int f[N];
f[0] = 1;
for (int i = 1; i <= n; ++i) {
f[i] = 0;
for (int j = 1; j <= i; ++j)
f[i] = add(f[i], 1ll * C[i - 1][j - 1] * f[j - 1] % Mod * f[i - j] % Mod);
f[i] = 1ll * f[i] * mi(x, i) % Mod;
}
return f[n];
}
inline void Lagrange() {
g[0] = 1;
for (int i = 1; i <= m; ++i)
for (int j = i; ~j; --j)
g[j] = add(1ll * (Mod - i) * g[j] % Mod, g[j - 1]);
for (int i = 1; i <= m; ++i) {
int mul = 1;
for (int j = 1; j <= m; ++j)
if (i != j)
mul = 1ll * mul * dec(i, j) % Mod;
mul = 1ll * y[i] * mi(mul, Mod - 2) % Mod;
for (int j = m - 1, res = g[m]; ~j; --j)
f[j] = add(f[j], 1ll * mul * res % Mod), res = add(g[j], 1ll * i * res % Mod);
}
}
signed main() {
scanf("%d%d", &n, &lim), m = n * (n + 1) / 2 + 1;
if (lim >= m)
return puts("0"), 0;
else if (n == 1)
return puts("1"), 0;
C[0][0] = 1;
for (int i = 1; i <= n; ++i) {
C[i][0] = 1;
for (int j = 1; j <= i; ++j)
C[i][j] = add(C[i - 1][j], C[i - 1][j - 1]);
}
for (int i = 1; i <= m; ++i)
y[i] = calc(i);
Lagrange();
int ans = 0;
for (int i = lim; i <= m; ++i)
ans = add(ans, f[i]);
printf("%d", ans);
return 0;
}
边双连通分量
给出一棵 \(n\) 个点的树，可以加入若干条边（不能形成重边和自环），对于每个 \(1 \leq i \leq n\) ，求加边之后恰有 \(i\) 个边双的方案数 \(\bmod (10^9 + 7)\) 。
\(n \leq 700\)
恰好的限制是困难的，考虑容斥，则需要求出钦定有若干条割边的方案数。
设 \(f_{u, i, j}\) 表示考虑 \(u\) 的子树，最上面的连通块大小为 \(i\) ，钦定了 \(j\) 条割边的方案数，考虑 \(f_{u, i, j}\) 与 \(f_{v, p, q}\) 的合并：
\[f_{u, i, j} \times f_{v, p, q} \to f'_{u, i + p, j + q} \\
f_{u, i, j} \times f_{v, p, q} \to f'_{u, i, j + q + 1}
\]
直接做是 \(O(n^4)\) 的，考虑插值法，令：
\[F_{u, i}(x) = \sum_{j = 0} F_{u, i, j} x^j
\]
则转移为：
\[F_{u, i}(x) \times F_{v, p}(x) \to F'_{i, j + p}(x) \\
F_{u, i}(x) \times F_{v, p}(x) \times x \to F'_{i, j}
\]
由于 \(\deg F = n - 1\) ，于是带入 \(n\) 个点值，做 \(n\) 次 \(O(n^2)\) 的 DP 求出点值之后 \(O(n^2)\) 反推系数即可，时间复杂度 \(O(n^3)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 7e2 + 7;
struct Graph {
vector<int> e[N];
inline void insert(int u, int v) {
e[u].emplace_back(v);
}
} G;
int C[N][N], pw[N * N], f[N][N], siz[N], g[N], y[N], p[N], h[N];
int n, m;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline int sgn(int n) {
return n & 1 ? Mod - 1 : 1;
}
inline void prework(int n) {
pw[0] = 1;
for (int i = 1; i <= n * n; ++i)
pw[i] = 2ll * pw[i - 1] % Mod;
C[0][0] = 1;
for (int i = 1; i <= n; ++i) {
C[i][0] = 1;
for (int j = 1; j <= i; ++j)
C[i][j] = add(C[i - 1][j], C[i - 1][j - 1]);
}
}
void dfs(int u, int fa, int x0) {
siz[u] = 1, f[u][1] = 1;
for (int v : G.e[u]) {
if (v == fa)
continue;
dfs(v, u, x0), memset(g, 0, sizeof(int) * (siz[u] + siz[v] + 1));
for (int i = 0; i <= siz[u]; ++i)
for (int j = 0; j <= siz[v]; ++j)
g[i + j] = add(g[i + j], 1ll * f[u][i] * f[v][j] % Mod);
memcpy(f[u], g, sizeof(int) * (siz[u] + siz[v] + 1)), siz[u] += siz[v];
}
for (int i = 1; i <= siz[u]; ++i)
f[u][0] = add(f[u][0], 1ll * f[u][i] * pw[i * (i - 1) / 2 - (i - 1)] % Mod * x0 % Mod);
}
namespace Lagrange {
int f[N], g[N];
inline void solve(int *y, int n) {
g[0] = 1;
for (int i = 1; i <= n; ++i)
for (int j = i; ~j; --j)
g[j] = add(1ll * (Mod - i) * g[j] % Mod, g[j - 1]);
for (int i = 1; i <= n; ++i) {
int mul = 1;
for (int j = 1; j <= n; ++j)
if (i != j)
mul = 1ll * mul * dec(i, j) % Mod;
mul = 1ll * y[i] * mi(mul, Mod - 2) % Mod;
for (int j = m - 1, res = g[n]; ~j; --j)
f[j] = add(f[j], 1ll * mul * res % Mod), res = add(g[j], 1ll * i * res % Mod);
}
}
} // namespace Lagrange
signed main() {
scanf("%d", &n), m = n + 1;
for (int i = 1; i < n; ++i) {
int u, v;
scanf("%d%d", &u, &v);
G.insert(u, v), G.insert(v, u);
}
prework(n);
for (int i = 1; i <= m; ++i)
memset(f, 0, sizeof(f)), dfs(1, 0, i), y[i] = f[1][0];
Lagrange::solve(y, m);
for (int i = 1; i <= n; ++i) {
int ans = 0;
for (int j = i; j <= n; ++j)
ans = add(ans, 1ll * sgn(j - i) * C[j - 1][i - 1] % Mod * Lagrange::f[j] % Mod);
printf("%d ", ans);
}
return 0;
}
浙公网安备 33010602011771号