拉格朗日反演

What · 是什么

拉格朗日反演（Lagrange Inversion）给出形式幂级数复合逆系数的闭式：若 f(x)=Σ_{k>=1} f_k x^k 且 f_1≠0，则 f 有唯一复合逆 g（即 f(g(x))=x）。对任意形式幂级数（或形式 Laurent 级数）Φ 有 [x^n]Φ(f(x))=[x^{n-1}]Φ(x)·g'(x)/g(x)·(x/g(x))^n。特例 [x^n] f(x)^k = (k/n)[x^{n-k}](x/g(x))^n。利用该公式可在 O(n log n)（NTT 加速幂与导数）内求复合逆的前 n 项。是生成函数计数（树、平面二叉树、带标号组合结构等）的高级工具，用于求解'自指方程'的显式系数。

Why · 为什么

How · 实现

伪代码

# 求 f 的复合逆 g 前 n 项
# 用 [x^n]g(x) = (1/n)[x^{n-1}](x/f(x))^n
inv_f = inverse(f/x, n)      # 幂级数逆
for n in 1..N:
    g[n] = (1/n) * [x^{n-1}] (inv_f)^n   # 幂用 NTT 快速幂
# 或者逐项解 f(g)=x(朴素)

Python 实现

MOD = 998244353


def lagrange_inversion_naive(f, n):
    """拉格朗日反演(朴素 O(n^3), 教学验证版):
    f(x) = sum_{k>=1} f_k x^k, f_1 != 0. 求复合逆 g 使 f(g(x)) = x, 前 n 项.
    公式: [x^n] g(x) = 1/n * [x^{n-1}] (x/f(x))^n.
    这里用幂级数除法直接解 g: f(g(x)) = x 逐项确定 g 系数。"""
    f = f[:n]
    assert f[0] == 0 and f[1] != 0
    g = [0] * n
    g[1] = pow(f[1], MOD - 2, MOD)
    # f(g(x)) = sum f_k g(x)^k = x
    # 逐项: [x^n] sum_{k>=1} f_k g(x)^k = delta_{n,1}
    for m in range(2, n):
        # 未知 g[m], 贡献来自 k=1 项: f_1 * g[m] * x^m
        # 已知项: sum_{k>=2} f_k [x^m] g(x)^k + f_1*[x^m] g(x) 中 g_m 未知
        # 用递推: f_1 g_m + sum_{k>=2} f_k [x^m] g^k = 0
        s = 0
        for k in range(2, n):
            # [x^m] g^k, 只含 g_0..g_{m-1} (因为 k>=2, g_1>=1 次)
            c = coeff_pow(g, k, m)
            s = (s + (f[k] if k < len(f) else 0) * c) % MOD
        g[m] = (-s * pow(f[1], MOD - 2, MOD)) % MOD
    return g


def coeff_pow(g, k, m):
    """[x^m] g(x)^k, g[0]=0, 用卷积逐层。O(k*m^2) 教学版。"""
    cur = [1]
    for _ in range(k):
        new = [0] * (len(cur) + len(g) - 1)
        for i, x in enumerate(cur):
            if x == 0:
                continue
            for j, y in enumerate(g):
                new[i + j] = (new[i + j] + x * y) % MOD
        cur = new
        if len(cur) > m + 1:
            cur = cur[:m + 1]
    return cur[m] if m < len(cur) else 0


def compose_verify(g, f, n):
    """验证 f(g(x)) = x mod x^n。"""
    # g 已是 f 的复合逆, 验证 f(g) 用逐项合成
    cur = [1]
    res = [0] * n
    for k in range(n):
        if k == 0:
            cur = [1]
        else:
            new = [0] * (len(cur) + len(g) - 1)
            for i, x in enumerate(cur):
                if x == 0:
                    continue
                for j, y in enumerate(g):
                    if i + j < n:
                        new[i + j] = (new[i + j] + x * y) % MOD
            cur = new[:n]
        c = f[k] if k < len(f) else 0
        if c:
            for i in range(min(len(cur), n)):
                res[i] = (res[i] + c * cur[i]) % MOD
    return res


if __name__ == "__main__":
    # 例: f(x) = x + x^2, 复合逆满足 g + g^2 = x
    # g(x) = x - x^2 + 2x^3 - 5x^4 + ... (Catalan 数相关)
    f = [0, 1, 1]
    g = lagrange_inversion_naive(f, 10)
    # 验证 f(g(x)) = x
    c = compose_verify(g, f, 10)
    assert c == [0, 1, 0, 0, 0, 0, 0, 0, 0, 0], c
    assert g[0] == 0 and g[1] == 1 and g[2] == MOD - 1, g  # x - x^2 + ...
    # 例: f(x) = x - x^2/2, 复合逆 g + g^2/2 = x
    inv2 = (MOD + 1) // 2
    f2 = [0, 1, (-inv2) % MOD]
    g2 = lagrange_inversion_naive(f2, 8)
    c2 = compose_verify(g2, f2, 8)
    assert c2 == [0, 1, 0, 0, 0, 0, 0, 0], c2
    print("lagrange-inversion tests passed")

易错点

• f(0)≠0 或 f'(0)=0 时无复合逆还硬算
• 留数/链式换元公式中 ord(g) 因子遗漏
• [x^n]f^k 公式中的 k/n 系数漏乘
• 形式 Laurent 级数的负次项处理错误
• 朴素逐项法只到小规模, 大规模需 NTT 库

例题

信源

• OI Wiki✓ https://oi-wiki.org/math/poly/lagrange-inversion/
• OI Wiki✓ https://raw.githubusercontent.com/OI-wiki/OI-wiki/master/docs/math/poly/lagrange-inversion.md