拉格朗日插值

拉格朗日插值

一般性插值

P4781 【模板】拉格朗日插值

给出 \(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;
}
posted @ 2025-07-20 10:53  wshcl  阅读(17)  评论(0)    收藏  举报