多项式牛顿迭代

一、引言：推开高效多项式算法的大门

在算法竞赛的数学领域，多项式占据着举足轻重的地位。无论是快速傅里叶变换 (FFT) / 数论变换 (NTT) 在加速卷积运算中的惊艳表现，还是生成函数在组合计数问题中的巧妙运用，都离不开对多项式各种运算的深刻理解和熟练掌握。当我们面对一些更复杂的多项式方程，例如求解多项式的逆元、平方根、对数、指数，乃至解决一些特定形式的函数方程时，一个强大而普适的工具——多项式牛顿迭代——便应运而生。

你可能对数值分析中的牛顿迭代法有所耳闻，它是一种求解非线性方程 根的有效迭代方法。其核心思想是利用函数在某一点的泰勒展开，通过线性逼近来逐步靠近方程的根。多项式牛顿迭代法正是这一思想在多项式域上的自然推广。它的目标是求解形如 的方程，其中 是一个关于多项式 的“函数”，我们希望找到一个多项式 满足这个同余方程。

掌握多项式牛顿迭代，意味着你将能够优雅地解决一系列看似棘手的多项式问题。

二、温故知新：数值牛顿迭代法

在深入多项式牛顿迭代之前，让我们先回顾一下它在数值分析中的原型——牛顿-拉弗森方法（Newton-Raphson method），通常简称为牛顿迭代法。

假设我们需要求解非线性方程 的一个根。牛顿迭代法的步骤如下：

1. 选取一个初始猜测值 。

2. 进行迭代：

   其中 是函数 在 处的导数。

这个迭代公式是如何得到的呢？考虑函数 在点 处的一阶泰勒展开：

我们希望找到一个 使得 ，所以令上式右边为0：

解出 ，得到：

我们将这个新解出的 作为下一次迭代的近似值 。

几何上，牛顿迭代法可以解释为：在点 处作函数 的切线，该切线与 轴的交点即为 。

牛顿迭代法的一个显著优点是其收敛速度。如果初始猜测值 足够接近真实根，并且 在根附近不为零，那么牛顿迭代法具有二次收敛性，即每次迭代后，有效数字的位数大约会翻倍。这意味着它能非常迅速地逼近方程的根。

这个核心思想——通过泰勒展开进行线性逼近，并迭代求解——将被我们巧妙地移植到多项式的世界。

三、核心思想：多项式域上的牛顿迭代

现在，我们将牛顿迭代的思想推广到多项式环上。我们要解决的问题通常是这样的：给定一个关于多项式 的函数方程 ，求解多项式 。这里的 是一个“作用于多项式”的算子，例如，（求逆元 使得 ），或者 （求平方根 使得 ）。

假设我们已经求得一个多项式 ，它在模 意义下满足方程，即：

我们的目标是找到一个多项式 ，使得它在模 (或者某个更高的幂次，通常是 ) 意义下满足方程：

并且 。这意味着 和 的前 项系数是相同的。

为了从 得到 ，我们对 在 处进行泰勒展开。注意，这里的泰勒展开是针对多项式 的，而不是针对 的。我们将 看作是一个以多项式 为自变量的”函数”。其泰勒展开（只取前两项）可以形式化地写作：

这里的 表示 对其多项式参数 在 处的“形式导数”或“函数导数”。这个概念初看起来可能有些抽象，我们会在具体应用中看到它的实际含义。

我们希望 ，所以：

从这个式子中解出 :

于是得到多项式牛顿迭代的公式：

这个公式与数值牛顿迭代的公式 在形式上是完全一致的！

关键性质与分析：

1. 倍增性质：如果 ，那么 的最低次项至少是 。
   同时，我们有 ，这意味着 的系数从 到 都为0。
   因此， 中的项实际上是从 到 的。
   当 是一个常数或者其逆元存在且良好定义时，整个迭代过程是有效的。

2. 收敛性：迭代的次数是对数级别的。如果我们要计算到 ，初始时我们通常有一个 (即常数项) 的解 (即 )。然后通过迭代得到 , , , …, 直到 其中 。这个过程需要 次迭代。

3. 复杂度：在每次迭代中，计算 和 可能涉及到多项式乘法、加法等操作。如果当前处理的长度为 ，那么多项式乘法（使用NTT）的复杂度是 。由于 每次翻倍，总的复杂度是 。这是一个收敛的级数，总复杂度为 ，其中 是最终需要的项数。

迭代的起点 (Base Case):
通常，我们需要一个初始的 使得 。这通常意味着求解 的常数项 。令 ，则 。这个 往往可以直接从原方程中解出。例如，若求 的逆元 ，则 。令 ，有 ，所以 (模意义下)。

四、多项式运算的基石：NTT与基本操作

在深入探讨牛顿迭代的具体应用之前，我们需要先准备好一些基本的多项式运算工具。这些运算本身就是多项式算法的基石，并且它们的实现效率直接影响牛顿迭代的总效率。我们将使用数论变换 (NTT) 来实现多项式乘法。

4.1 数论变换 (NTT)

NTT 是 FFT 在有限域上的变体，它利用模意义下的原根替代了复数单位根，从而避免了浮点精度问题，非常适合在高水平程序设计竞赛中使用。

#include<bits/stdc++.h>

using namespace std;
using ll = long long;

const int MOD = 998244353; // 常用NTT模数
const int G = 3;           // MOD的原根
const int IG = (MOD + 1) / G; // G在模MOD下的逆元

// 快速幂
ll qpow(ll b, ll e) {
    ll r = 1;
    b %= MOD;
    while (e > 0) {
        if (e & 1) r = r * b % MOD;
        b = b * b % MOD;
        e >>= 1;
    }
    return r;
}

ll inv(ll n) {
    return qpow(n, MOD - 2);
}

using poly = vector<ll>; // 定义多项式类型

// NTT预处理变址
vector<int> g_rev; // 全局或传递的rev数组
void pre_ntt(int lim) {
    if ((int)g_rev.size() == lim) return; // 避免重复计算
    g_rev.assign(lim, 0);
    for (int i = 0; i < lim; ++i) {
        g_rev[i] = (g_rev[i >> 1] >> 1) | ((i & 1) ? (lim >> 1) : 0);
    }
}

// NTT核心实现
// op = 1 为NTT, op = -1 为INTT
void ntt(poly& p, int lim, int op) {
    if ((int)p.size() != lim) p.resize(lim); // 补全或调整长度
    for (int i = 0; i < lim; ++i) {
        if (i < g_rev[i]) swap(p[i], p[g_rev[i]]);
    }

    for (int mid = 1; mid < lim; mid <<= 1) {
        ll wn = qpow(op == 1 ? G : IG, (MOD - 1) / (mid << 1));
        for (int j = 0; j < lim; j += (mid << 1)) {
            ll w = 1;
            for (int k = 0; k < mid; ++k, w = w * wn % MOD) {
                ll x = p[j + k];
                ll y = w * p[j + mid + k] % MOD;
                p[j + k] = (x + y) % MOD;
                p[j + mid + k] = (x - y + MOD) % MOD;
            }
        }
    }

    if (op == -1) {
        ll ilim = inv(lim);
        for (int i = 0; i < lim; ++i) {
            p[i] = p[i] * ilim % MOD;
        }
    }
}

// 多项式乘法
poly mul(poly p1, poly p2) {
    int s1 = p1.size(), s2 = p2.size();
    if (s1 == 0 || s2 == 0) return {};
    int lim = 1;
    while (lim < s1 + s2 - 1) lim <<= 1;

    pre_ntt(lim); // 预处理变址

    ntt(p1, lim, 1);
    ntt(p2, lim, 1);

    poly res(lim);
    for (int i = 0; i < lim; ++i) {
        res[i] = p1[i] * p2[i] % MOD;
    }

    ntt(res, lim, -1);
    res.resize(s1 + s2 - 1); // 调整到实际长度
    return res;
}

• 时间复杂度: ntt 为 ，mul 为 ，其中 是扩展到的2的幂次长度。
• 空间复杂度: 。 g_rev 也是 。

4.2 多项式加减与数乘

这些操作比较简单，直接对系数进行操作即可。

// 多项式加法
poly add(poly p1, poly p2) {
    int s = max(p1.size(), p2.size());
    p1.resize(s); p2.resize(s); // 对齐长度
    poly r(s);
    for (int i = 0; i < s; ++i) {
        r[i] = (p1[i] + p2[i]) % MOD;
    }
    return r;
}

// 多项式减法
poly sub(poly p1, poly p2) {
    int s = max(p1.size(), p2.size());
    p1.resize(s); p2.resize(s); // 对齐长度
    poly r(s);
    for (int i = 0; i < s; ++i) {
        r[i] = (p1[i] - p2[i] + MOD) % MOD;
    }
    return r;
}

// 多项式数乘
poly smul(poly p, ll k) {
    k = (k % MOD + MOD) % MOD; // 保证k为[0, MOD-1]
    poly r = p;
    for (ll &x : r) {
        x = x * k % MOD;
    }
    return r;
}

// 多项式截断 (取前deg项)
void redeg(poly& p, int deg) {
    if (deg < 0) deg = 0;
    if ((int)p.size() > deg) {
        p.resize(deg);
    } else if ((int)p.size() < deg) { // 可选：如果需要补0到deg项
        p.resize(deg, 0);
    }
}

4.3 多项式求逆 (Polynomial Inverse)

给定多项式 , 求 使得 。
迭代公式: 。
其中 是 的解。

Base Case: 当 (即 ) 时, 。所以 。这要求 。

// 题意概括: 计算多项式 P(x) 在模 x^deg 意义下的逆元。
// P(x) * Q(x) = 1 (mod x^deg)
poly pInv(poly p, int deg) {
    if (deg == 0) return {};
    poly q = {inv(p[0])}; // 基准情况: q_0 = p_0^{-1}
    int kl = 1; // 当前长度 k

    while (kl < deg) {
        kl <<= 1;
        poly pa = p; 
        redeg(pa, kl); // P(x) mod x^{kl}

        int lim = 1;
        while (lim < kl + kl) lim <<= 1; // NTT长度，A*B0^2需要2*kl，这里A*B0，B0*(...)
                                         // A(kl项), B0(kl/2项) -> A*B0 (kl+kl/2-1项)
                                         // B0(kl/2项), (2-AB0)(kl项) -> B0*(2-AB0) (kl+kl/2-1项)
                                         // 所以lim至少是kl+kl/2，或者更宽松取2*kl (2*kl覆盖所有中间步骤)
        pre_ntt(lim);

        poly npa = pa;  // pa的NTT形式
        poly nq = q;    // q (即B0)的NTT形式
        ntt(npa, lim, 1);
        ntt(nq, lim, 1);

        for (int i = 0; i < lim; ++i) {
            // nq = nq * (2 - npa * nq) % MOD
            ll vq = nq[i];
            ll vpaq = npa[i] * vq % MOD; // A * B0 (点值)
            nq[i] = vq * (2 - vpaq + MOD) % MOD; // B0 * (2 - A*B0) (点值)
        }
        
        ntt(nq, lim, -1); // 逆NTT回系数
        q = nq;
        redeg(q, kl); // 截断到 x^{kl}
    }
    redeg(q, deg); // 最终截断到 x^deg
    return q;
}
// 复杂度: T(N) = T(N/2) + O(N log N) => O(N log N)
// 空间复杂度: O(N)

4.4 多项式微分与积分 (Polynomial Derivative and Integral)

对于多项式 :
其导数 。
其积分 。积分常数 通常为0。

// 题意概括: 计算多项式 P(x) 的导数。
poly pDer(poly p) {
    if (p.empty()) return {};
    int s = p.size();
    if (s == 1) return {0}; // 常数导数为0
    poly r(s - 1);
    for (int i = 1; i < s; ++i) {
        r[i-1] = (ll)i * p[i] % MOD;
    }
    // 如果原多项式是0，结果是空vector (pDer({0}) -> {})
    // 这可以通过约定处理，例如pDer({0}) -> {0}。当前实现会是{}
    // 如果pDer({c}) c!=0 -> {0}
    if (r.empty() && !p.empty() && p[0] != 0 && s == 1) return {0}; // 明确处理常数
    if (r.empty() && s > 1) return {0}; // 例如 P(x) = x，P'(x) = 1，而不是空
    return r;
}
// 复杂度: O(S)，S为多项式项数
// 空间复杂度: O(S)

// 预计算模逆元 (用于积分)
vector<ll> itbl; // inv_tbl
void pre_itbl(int s) {
    if (itbl.empty()) itbl.push_back(0); // itbl[0] 无意义
    int os = itbl.size();
    if (s >= os) {
        itbl.resize(s + 1);
        for (int i = os; i <= s; ++i) {
            if (i == 0) continue; // 1/0无意义
            itbl[i] = inv(i);
        }
    }
}

// 题意概括: 计算多项式 P(x) 的积分，常数项为0。
poly pInt(poly p) {
    if (p.empty()) return {0}; // 0的积分是C, 此处为0
    int s = p.size();
    pre_itbl(s + 1); // 需要 inv_tbl[1] 到 inv_tbl[s]
    poly r(s + 1);
    r[0] = 0; // 积分常数项为0
    for (int i = 0; i < s; ++i) {
        if (i + 1 > s) pre_itbl(i + 1); // 确保itbl够大
        r[i+1] = p[i] * itbl[i+1] % MOD;
    }
    return r;
}
// 复杂度: O(S) (预处理逆元后)
// 空间复杂度: O(S)

五、牛顿迭代的“小试牛刀”：多项式开方、对数与指数

有了前面的基础，我们现在可以运用牛顿迭代法来解决一些更复杂的多项式运算了。

5.1 多项式平方根 (Polynomial Square Root)

给定多项式 ，求 使得 。
令 。则 。
牛顿迭代公式: 。

Base Case: 。通常 , 则 。若 不是1，需解二次剩余。这里假设 。

// 题意概括: 计算多项式 A(x) 在模 x^deg 意义下的平方根 B(x)。
// B(x)^2 = A(x) (mod x^deg)。通常要求 A[0] = 1, 则 B[0] = 1。
poly pSqrt(poly pa, int deg) {
    if (deg == 0) return {};
    // 此处简化假设 pa[0]=1, 所以 b[0]=1
    // 若 pa[0] != 1, 需要二次剩余算法求 b[0]
    // if (pa.empty() || pa[0] != 1) { /* 错误或复杂情况处理 */ }
    
    poly pb = {1}; // B_0 = 1 (因为 A[0]=1)
    int kl = 1;    // 当前长度 k
    ll inv2 = inv(2);

    while (kl < deg) {
        kl <<= 1;
        poly ca = pa; // 当前A(x)
        redeg(ca, kl); // A(x) mod x^{kl}

        // B_new = (B_old + A / B_old) / 2 mod x^{kl}
        // A / B_old = A * inv(B_old)
        // inv(B_old) 需要是 B_old 在 mod x^{kl} 意义下的逆
        poly inv_pb = pInv(pb, kl); // pb是(kl/2)项, 求其 mod x^{kl} 的逆
        
        poly tdb = mul(ca, inv_pb); // A * B_0^{-1}
        redeg(tdb, kl);             // mod x^{kl}

        pb.resize(kl, 0); // 扩展pb到kl长度
        
        pb = add(pb, tdb);    // B_0 + A * B_0^{-1}
        pb = smul(pb, inv2);  // (B_0 + A * B_0^{-1}) / 2
        redeg(pb, kl);        // 确保长度正确
    }
    redeg(pb, deg);
    return pb;
}
// 复杂度: T(N) = T(N/2) + O(N log N) (主要来自pInv和mul) => O(N log N)
// 空间复杂度: O(N)

5.2 多项式对数 (Polynomial Logarithm)

给定多项式 且 ，求 。
利用 。
所以 。
然后 。
积分常数 由 确定。

// 题意概括: 计算多项式 A(x) 的对数 ln(A(x)) (mod x^deg)。
// 要求 A[0] = 1。
poly pLn(poly pa, int deg) {
    if (deg == 0) return {};
    // B'(x) = A'(x) * inv(A(x))
    // A'(x) 有 deg-1 项 (若A有deg项). inv(A(x)) 也需相应项数.
    // B(x) 要 deg 项, B'(x) 要 deg-1 项.
    // A'(x) mod x^{deg-1}, A(x)^{-1} mod x^{deg-1}

    poly dpa = pDer(pa);     // A'(x)
    // inv(A(x)) mod x^{deg}. 如果A有deg项, A'有deg-1项. B'有deg-1项.
    // A^{-1} 应该取到 mod x^{deg-1} (或deg, 然后截断).
    poly ipa = pInv(pa, deg); 
                                 
    poly bp = mul(dpa, ipa); // B'(x) = A'(x) * A(x)^{-1}
    redeg(bp, deg - 1 > 0 ? deg -1 : 0); // B'(x) mod x^{deg-1} (即前 deg-1 项)
                                         // deg-1可能为0或负，确保非负

    poly pb = pInt(bp);      // B(x) = int B'(x) dx
    redeg(pb, deg);          // 确保是 deg 项
    return pb;
}
// 复杂度: O(N log N) (主要来自pInv和mul)
// 空间复杂度: O(N)

关于 pLn 中 pInv 的度数修正：
pDer(pa) 的结果长度最多为 pa.size()-1。pInv(pa, deg) 的结果长度为 deg。
mul(dpa, ipa) 之后，bp 会被截断到 deg-1。
这意味着 dpa 和 ipa 自身在乘法前，如果项数超过 deg-1 (近似)，那些高次项对最终结果 bp (截断后)没有贡献。
更精确地，dpa 需要到 deg-1 项，ipa (即 ) 也需要到 deg-1 项。
所以 pInv(pa, deg-1>0 ? deg-1:1) 传递给 pInv 是更紧凑的，但传递 deg 也是安全的，因为后续会截断。 polyInv(a, deg) 中 deg 指的是结果的项数。pInv(pa, deg) 然后在 mul 后 redeg(bp, deg-1)。这没有问题。

5.3 多项式指数 (Polynomial Exponential / Exp)

给定多项式 且 , 求 。
即 。
令 。则 。
牛顿迭代公式: 。

Base Case: 。所以 。

// 题意概括: 计算多项式 exp(A(x)) (mod x^deg)。
// 要求 A[0] = 0。则结果 B[0] = 1。
poly pExp(poly pa, int deg) {
    if (deg == 0) return {};
    // 基准情况: B_0 = 1 (因 A[0]=0)
    poly pb = {1}; 
    int kl = 1; // 当前长度 k

    while (kl < deg) {
        kl <<= 1;
        
        // B_new = B_old * (1 - ln(B_old) + A) (mod x^{kl})
        // A_sub = A mod x^{kl}
        // ln(B_old) (mod x^{kl}). B_old此时有 kl/2 项.
        poly lnb = pLn(pb, kl); // pLn(B_old (kl/2项), kl) -> ln(B_old) (kl项)
        
        poly ca = pa; 
        redeg(ca, kl); // A(x) mod x^{kl}

        // term = 1 - ln(B_old) + A (mod x^{kl})
        poly term = sub(ca, lnb); // A - ln(B_old)
        // term[0] += 1
        if (term.empty()) term.push_back(0); // 若A-lnB是0，确保有位置给常数项
        term[0] = (term[0] + 1 + MOD) % MOD; // 1 + (A - ln(B_old))
        
        // 确保term有kl项，高位补0 (若A或lnB项数不足kl，sub后可能不足)
        // redeg在之前已经处理了ca和lnb的长度到kl (pLn返回kl项)
        // sub会按最长处理，所以term应该已经是kl项。
        // 如果ca原始长度小于kl，redeg(ca,kl)会补0。
        // 所以term应该有kl项。以防万一，可以redeg(term, kl)。
        redeg(term, kl);

        // B_new = B_old * term
        // pb (B_old) 有 kl/2 项. term 有 kl 项.
        pb = mul(pb, term); 
        redeg(pb, kl);
    }
    redeg(pb, deg);
    return pb;
}
// 复杂度: T(N) = T(N/2) + O(N log N) (pLn内部有pInv) => O(N log N)
// 空间复杂度: O(N)

这些基本运算——求逆、开方、对数、指数——构成了解决更复杂问题的“代数积木”。掌握它们的推导和高效实现至关重要。尤其是 pInv 和 pExp (其内部依赖 pLn，pLn 内部依赖 pInv)，它们是许多生成函数计数问题的核心。

六、牛顿迭代与生成函数：组合计数的利器

生成函数是组合数学中分析序列、解决计数问题的强大工具。一个序列 ${a_n}{n \ge 0}A(x) = \sum{n=0}^\infty a_n x^nA(x) = \sum_{n=0}^\infty \frac{a_n}{n!} x^n$。许多组合结构的计数问题可以转化为对其生成函数满足的某个方程的求解。多项式牛顿迭代法为此类方程的求解提供了有效的途径。

6.1 从递推关系到函数方程

许多组合对象的计数满足递推关系，这些递推关系往往可以直接翻译成其生成函数的代数方程或微分方程。

例如，考虑Catalan数的一个定义：, for 。
令 为Catalan数的OGF。 是 中 的系数。
所以 。
整理得 。
这是一个关于 的二次方程。我们可以用二次求根公式解出 。根据 , 我们选择合适的符号（通常是负号）并处理 的情况得到 。
这个求解过程涉及到多项式平方根。

更一般地，如果一个组合结构的生成函数 满足某个形式的函数方程 ，我们就可以尝试用牛顿迭代来求解 。

6.2 用牛顿迭代求解

假设我们要解方程 中的 ，其中 是一个关于多项式 和 的表达式。我们将 看作是 的函数，而 及其多项式 (如 ) 被视为系数。
牛顿迭代公式为：

其中 是模 的解， 是 对 的形式偏导数。

示例：求解
这里 是一个给定的多项式 (例如，在Catalan数的例子中 ）。
令 。
则 。
迭代公式为：

Base Case: 。由 得到 。解出 。
例如，对于Catalan数 , 。所以 , 或 。
由于 , 。所以 。
那么 的常数项是 ，其逆元存在。

// 题意概括: 求解方程 B(x)^2 - B(x) + A(x) = 0 (mod x^deg)
// B[0] 由 A[0] 决定: B[0]^2 - B[0] + A[0] = 0.
// 此处假设 A[0]=0, B[0]=0 或 B[0]=1. 根据具体问题选择.
// 例如卡特兰数 C(x) = xC(x)^2 + 1 -> xC^2 - C + 1 = 0, A(x)=x. C[0]=1.
// 本例用 T(x)^2 - T(x) + x = 0, T[0]=0. (T(x) = x + x^2 + 2x^3 + ...)
poly solveQuadFE(poly pa, int deg, ll b0) {
    if (deg == 0) return {};
    poly pb = {b0}; // B_0 = b_0
    int kl = 1;

    while (kl < deg) {
        kl <<= 1;
        
        poly ca = pa; // A(x)
        redeg(ca, kl); // A(x) mod x^kl
        
        // Numerator: pb^2 - pb + ca
        poly pb2 = mul(pb, pb);
        redeg(pb2, kl);
        
        poly num = sub(pb2, pb);
        num = add(num, ca);
        redeg(num, kl);

        // Denominator: 2*pb - 1
        poly den = smul(pb, 2);
        if (den.empty()) den.push_back(0); // 2*0=0
        den[0] = (den[0] - 1 + MOD) % MOD;
        redeg(den, kl); // Ensure length for inverse

        poly inv_den = pInv(den, kl); // (2*pb - 1)^{-1}

        poly term = mul(num, inv_den);
        redeg(term, kl);

        pb = sub(pb, term); // pb - term
        redeg(pb, kl);
    }
    redeg(pb, deg);
    return pb;
}
// 复杂度: T(N) = T(N/2) + O(N log N) => O(N log N)
// 空间复杂度: O(N)

在 solveQuadFE 中，如果 b0 使得 ，则分母的逆元不存在，迭代失败。这通常意味着方程的解在该点有奇点或不满足牛顿迭代的前提。

七、多项式任意次幂 (Polynomial Power)

计算 ，其中 可以是整数、分数，甚至在某些情况下是另一个多项式（非常数）。最通用的方法是利用对数和指数函数：
。

处理步骤与条件：

1. :
   可以直接计算 ，然后 (如果 是常数，这是数乘；如果 是多项式 ，这是多项式乘法 )，最后 。所有运算都在 下进行。
   复杂度为两次核心操作（pLn, pExp）和一次乘法（如果 是多项式），即 。

2. 且 :
   令 。 .
   的常数项是 。令 。
   则 。
   可以用快速幂计算。
   可以用上述 的方法计算 。
   这一步需要 在模 下有意义的 次幂。例如，如果 ，则 必须是二次剩余。

3. :
   设 , 其中 。
   令 。则 。
   。
   可以用第2种情况处理：。
   这个结果 的各项系数需要左移 位。
   如果 , 那么 (除非 即 或 )。
   这要求 是整数，或者 是使得 为整数的分数。如果 是非负整数，这总是可行的。如果 是负整数，则 必须可逆（即 ）。

// 题意: 计算 pa(x)^k (mod x^deg). k为ll常数.
poly pPow(poly pa, ll k, int deg) { // pa: A(x)
    if (deg == 0) return {}; // 若所需度数为0, 返回空

    int pz = 0; // pa(x)中x的最低次项的指数, 即pa(x) = x^pz * ph(x)中的pz
    while (pz < (int)pa.size() && pa[pz] == 0) { // 找到第一个非零系数的索引
        pz++;
    }

    // 特殊情况: pa(x)是零多项式 (所有系数为0)
    if (pz == (int)pa.size()) { 
        if (k == 0) { // 0^0 通常约定为1
            poly r(deg, 0); if (deg > 0 && !r.empty()) r[0] = 1; return r;
        }
        // 对于 k > 0, 0^k = 0. 若 k < 0, 0的负数次幂无定义 (或无穷大).
        return poly(deg, 0); // 此处假设 k > 0 或题目保证合法性
    }
    
    // pa(x) = x^pz * ph(x), 其中 ph[0] != 0
    // pa(x)^k = x^{pz*k} * ph(x)^k

    // 检查 pa[0]=0 且 k<0 的情况 (即 pz > 0 且 k < 0)
    if (k < 0 && pz > 0) { // 出现 0 的负数次幂, 错误
        throw runtime_error("pPow: A[0]=0 (pz>0) and k<0");
    }

    // 计算最终结果 x^{pz*k} 的最低次项
    // pz*k 可能非常大, 使用 unsigned long long 避免溢出比较
    unsigned long long fmd = (unsigned long long)pz * k; // 最终结果中x的最低次项指数
    if (k > 0 && fmd >= (unsigned long long)deg) { // 若 x^{pz*k} 的次数已超出所需度数, 结果为0
        return poly(deg, 0);
    }

    poly ph; // ph(x) = pa(x) / x^pz, 即 pa(x) 右移 pz 位
    for (size_t i = pz; i < pa.size(); ++i) {
        ph.push_back(pa[i]);
    }
    
    int hDg = deg; // H(x)^k 需要计算到的项数 (模 x^hDg 意义下)
    if (k > 0) { // 若 k>0, 目标是 (x^pz H(x))^k mod x^deg.
                 // H(x)^k 自身需要计算到 mod x^{deg - pz*k} 
        if (fmd < (unsigned long long)deg) { // 确保 deg - pz*k 不为负
            hDg = deg - fmd;
        } else { // 结果应为0 (已被上面 fmd >= deg 的判断覆盖), 为hDg设安全值
            hDg = 0;
        }
    }
    // 若hDg为0, pLn/pExp等应能正确处理(通常返回空poly或含常数项的poly)
    redeg(ph, hDg); // 截断或扩展ph到hDg项

    // 处理ph为空但hDg>0的情况 (通常不应发生若pz和deg计算正确)
    if (ph.empty() && hDg > 0) { 
        // H(x)=0 但需要计算其幂的非0项.
        // 若 k=0, H(x)^0 = 1. 若 k>0, H(x)^k = 0.
        if (k==0) { poly r(deg,0); if(deg>0 && !r.empty()) r[0]=1; return r;}
        return poly(deg,0);
    }
    
    ll h0 = ph.empty() ? 1 : ph[0]; // H[0]. 若ph空(因hDg=0), h0取1不影响k=0时h0k=1
    if (ph.empty() && k > 0) h0 = 0; // 若ph空且k>0 (即H(x)=0), 则H[0]=0

    ll h0k; // H[0]^k
    if (k == 0) {
        h0k = 1; // 任何非零数(包括H[0])的0次幂是1. 0^0也约定为1.
    } else if (h0 == 0) { // H[0]=0 且 k!=0
        // 由前置检查, k此时必为正. 0^k = 0 for k>0.
        h0k = 0; 
    } else { // H[0] != 0 且 k != 0
        ll expB = k; // 实际指数基数
        if (k < 0) expB = -k; // 取绝对值

        ll effE = expB; // 有效指数 (模 MOD-1)
        if (MOD > 1 && (MOD-1) > 0) { // 仅当 MOD-1 > 0 时取模才有意义
           effE = expB % (MOD - 1);
           // 若 expB % (MOD-1) == 0 且 expB != 0 (即 k 不是0的倍数但整除MOD-1), 有效指数应为 MOD-1
           if (effE == 0 && expB != 0) effE = MOD - 1; 
        } else if (expB != 0 && MOD <= 1) { // MOD过小或无效时的特殊处理(竞赛不常见)
           // 此处简化，竞赛模数通常是大的素数
           // 对于 h0=1, 1^k = 1. h0=-1, (-1)^k. 其它 h0, 结果可能是0.
           if (h0 == 1) h0k = 1;
           else if (h0 == MOD -1 ) h0k = (effE % 2 == 0) ? 1 : (MOD -1); // h0 = -1
           else h0k = 0; // 其他情况在小MOD下简化为0
           // goto skip_qpow; // 跳过标准qpow
        }
        // 若 MOD-1 = 0 (即 MOD=1 或 MOD=2), expB % (MOD-1) 会出问题
        // MOD=2, MOD-1=1. expB % 1 = 0. 此时 effE = 1.
        if (MOD == 2 && expB !=0) effE = 1; // 在GF(2)中 a^e = a if a!=0,e>0

        if (k > 0) h0k = qpow(h0, effE);
        else h0k = qpow(inv(h0), effE); // k<0, H[0]^(-k) = (H[0]^{-1})^(-k)
    }
    // skip_qpow:; // 标签用于跳转

    poly rHk; // 结果 H(x)^k
    if (hDg == 0) { // H(x)^k mod x^0. 结果是常数h0k, 但多项式表示为0项 (或含1项{h0k}).
        rHk = {}; // 空多项式
    } else if (h0 == 0 && k > 0) { // H(x) 的常数项为0 (ph可能形如 {0, c1, ...}) 且 k>0
                                   // 这意味着 H(x) = x H'(x), H(x)^k = x^k H'(x)^k.
                                   // 但我们已分离 x^pz, 所以 ph[0] (即h0) 不应为0除非ph本身是0.
                                   // 此分支对应 ph = {0,0,...} (全零), 则 H(x)^k = 0.
        rHk.assign(hDg, 0);
    } else if (k == 0) { // H(x)^0 = 1 (只要H(x)非严格零多项式)
        rHk.assign(hDg, 0); if (hDg > 0) rHk[0] = 1;
    }
    else { // H[0]!=0, k!=0. 标准路径: H(x)^k = H[0]^k * exp(k * ln(H(x)/H[0]))
        poly phn = smul(ph, inv(h0)); // phn = ph / h0, 常数项为1
        redeg(phn, hDg);

        poly lnp = pLn(phn, hDg);   // lnp = ln(phn)
        // k*ln(phn). smul内部的k会模MOD, 这是对每个系数乘以k.
        // k可能很大或负数. (c_i * k) mod MOD.
        poly klp = smul(lnp, k);     
        redeg(klp, hDg); 
        // pExp 要求输入多项式常数项为0. ln(常数项1的多项式) 结果常数项为0.
        // k*0=0. 所以 klp[0] 应该是0.
        if(!klp.empty() && hDg > 0 && klp[0]!=0) {
            // 此情况理论上不应发生，若lnp[0]=0, 则k*lnp[0]=0
            // 若发生，强制设为0或检查逻辑
            klp[0]=0; 
        }

        poly expk = pExp(klp, hDg); // expk = exp(klp)
        
        rHk = smul(expk, h0k); // H(x)^k = h0k * expk
        redeg(rHk, hDg);
    }

    // 组合 x^{pz*k} * rHk 得到最终结果
    poly res(deg, 0); // 最终结果多项式
    if (fmd < (unsigned long long)deg) { // 如果x的最低次幂小于所需度数
        for (int i = 0; i < (int)rHk.size(); ++i) { // 将rHk的系数填充到res对应位置
            if (fmd + i < (unsigned long long)deg) { // 确保不越界
                res[fmd + i] = rHk[i];
            } else { // 超出所需度数，停止
                break;
            }
        }
    }
    return res;
}
// 复杂度: O(N log N), N=deg. 主要开销来自 pLn 和 pExp.
// 空间复杂度: O(N).

关于 pPow 中 k 的处理：
smul(poly P, ll val_k): 内部会将 val_k 对 MOD取模。
当计算 时，我们希望的是每个系数 乘以 。
即 。
所以 smul(ln_poly, k) 是正确的，smul 内部的 k % MOD 也是正确的，因为我们只关心系数模 MOD 的结果。
h0k = qpow(h0, k): 这里的 作为指数。如果 是素数，根据费马小定理， (若 )。所以指数 可以对 取模 (如果 )。如果 ，则 应该是 。qpow 通常能处理负指数（通过先求逆元）。
我们的 qpow 实现只处理正指数 e。若要支持负指数 ，应计算 qpow(inv(h0), -k)。
对于 qpow(base, exp)，exp 通常取模 MOD-1 (如果 MOD 是素数且 base != 0) 或 phi(MOD)。
在 pPow 的代码中，qpow(h0, k) 应该传入真实的 。如果 很大或负，qpow 需要正确处理。目前的 qpow 只接受 ll e >= 0。
为支持负 (且 ): qpow(h0, k) 应该改为 k>=0 ? qpow(h0, k) : qpow(inv(h0), -k)。
当 很大时，e 在 qpow 中是 k % (MOD-1) (或者对于h0=0的特殊处理)。
所以 h0k = qpow(h0, k) 这一行，如果 k 可能为负，需要调整 qpow 或调用方式。
在竞赛中，k 通常是非负整数。如果 可以是负数，意味着求 , 即 。
可以直接求 然后对其进行正整数次幂。

八、牛顿迭代解一般函数方程的推广

之前的多项式求逆、开方、指数函数都可以视为求解特定函数方程 的特例。牛顿迭代法的思想可以推广到更一般的函数方程，形如 ，只要我们能有效地计算 本身以及 对 的偏导数 $F’{B(x)}F’{B(x)}$ 的常数项可逆。

考虑方程 ，其中 是一个 “作用于多项式的函数” （例如 或 ）， 是给定多项式。
我们想求 。设 。
则迭代式为 。
这要求 及其导数 能够有效地作用于多项式。
例如，（一个固定的多项式P作用于Y）。
。这需要计算多项式的幂和总和。
。

一个更复杂的例子：
这在某些组合问题中出现，特别是涉及到映射或者特定类型的树的计数。
改写为 。
这里 被看作是 中的系数。
对 的形式导数是：

。

迭代公式为：

Base Case: 。由 。
如果 , 则 。此时 , 可逆。
如果 , 这个方程可能需要数值方法解 。在组合问题中，通常 的性质 (如 或 很小) 会使 容易确定。

复杂度分析：
在每次迭代中，我们需要计算 ，这本身是一个 的操作 (如果 是 项)。
然后是若干次多项式乘法和加减法。
总迭代次数是 次。如果每次迭代都是 ，那么总复杂度可能是 。
这是因为 从 变化， 。
不， 中的 是已知到 项的，计算 (或 ，取决于迭代项如何使用)。
计算 的复杂度是 。
在迭代 :

1. 计算 (或 就够了，取决于 和 如何依赖它)。
   ：这里 需到 ， 也需到 。
   本身是 项的。 是一个 项多项式。
   如果用 来计算所有项，那么 也是 的。
   则 (分子) 和 (分母) 都是 的。
   修正项 也是 。
   。这意味着 correction 的项必须准确到 。
   所以 和 都需要计算到 。
   (到 )， (到 )。
   : 这是 。
   然后若干 的乘法。
   所以每次迭代的复杂度是 其中 是当前的目标长度 ( in notation above)。
   总复杂度依然是 。

九、实战技巧与注意事项

应用牛顿迭代及其相关多项式算法时，一些细节和常见问题需要注意：

9.1 常数项处理

多项式开方、对数、指数等操作对其参数的常数项有特定要求：

• pInv(A, deg): 必须在模 下可逆 (即 )。
• pSqrt(A, deg): 必须是模 的二次剩余。竞赛中通常 , 则 (或 , 取1)。如果 , 则 。如果 , 需要先提取 的因子, , 。
• pLn(A, deg): 要求 。如果 , 。 是一个常数，但 可能没有良好定义（离散对数难求）。因此，通常题目会保证 或可以通过变换达到。
• pExp(A, deg): 要求 。如果 , 。 的常数项为0。 可以通过 计算，如果 足够小，或者 有其它方式计算。

9.2 度数控制 (redeg)

在迭代的每一步，多项式的长度（度数）需要精确控制。

• NTT 的长度 lim 必须是大于等于结果多项式度数加1的最小2的幂次。
• 例如，mul(p1, p2) 后结果长度是 p1.size() + p2.size() - 1。
• pInv, pSqrt, pExp 等函数在迭代时，当前处理的多项式 长度为 ，目标 长度为 。所有中间计算涉及的多项式（如 , ）都应截断或扩展到长度 进行运算，然后结果再截断到 。redeg(poly, length) 函数至关重要。不正确的长度处理是常见的bug来源。

9.3 模数选择与NTT的限制

• NTT 要求模数 是一个质数，且 有大的2的幂次因子，以支持足够长的变换。如 。
• 如果题目给定的模数不是NTT友好型，可能需要：
  • 中国剩余定理 (CRT)：选择多个NTT模数计算，然后合并结果。实现复杂，常数大。
  • 任意模数FFT/NTT：基于复数FFT或更复杂代数构造的技巧，通常更慢。
  • 观察题目数据范围是否允许 算法。

9.4 调试多项式代码

多项式相关代码（尤其是NTT和牛顿迭代）调试难度较大。

• 单元测试：对每个基本操作 (mul, pInv, pLn, pExp, pSqrt 等) 编写独立的测试用例。
  • 例如：。pInv({1, MOD-1}, 5) 应得到 {1,1,1,1,1}。
  • 。测试 pLn 时，输入 (常数项为1)，应得到 。
  • 。测试 pExp({0,1}, 5) (输入 ) 应得到 。
• 打印中间结果：在牛顿迭代的循环中，打印 和计算出的 的前几项，观察是否按预期收敛。
• 小数据模拟：手动模拟小规模数据下的迭代过程，与程序输出对比。
• 注意浮点精度问题（如果用FFT）与整型溢出。

9.5 复杂度陷阱

• 确保所有NTT长度 lim 的计算正确，避免不必要的过大长度。
• 牛顿迭代的递归式 保证了总复杂度 。要小心在 部分不慎引入了额外因子，比如对一个长度 的多项式进行了 次 操作，导致 。
• 例如，在 pExp 的迭代中，计算 ，如果 长 ，pLn(B_0, k) 的复杂度是 。这是正确的。

十、典型例题选讲

掌握理论后，通过题目实践是提升的关键。

例1：特定结构的树计数 (如Cayley’s Formula相关或更简单结构)

题意概括: 求解 ，其中 。这对应某些类型平面树的生成函数，其中 的系数是大小为 的树的数量。 (Catalan数生成函数 与此密切相关， 如果定义 为 那么 其中 是 的解)。

解法:
此方程 。
前面在6.2节中已经给出了这类方程的迭代求解框架 solveQuadFE。
这里 。 (因为 )。
所以调用 solveQuadFE({0,1}, N, 0) 即可。（{0,1} 代表多项式 ）。

#include<bits/stdc++.h>

using namespace std;
using ll = long long;

const int MOD = 998244353;
const int G = 3;
const int IG = (MOD + 1) / G;

// ... (qpow, inv, poly type, g_rev, pre_ntt, ntt) ...
// ... (add, sub, smul, redeg) ...
// ... (pInv, pDer, pre_itbl, pInt, pLn, pExp, pSqrt) ...
// 这些函数定义如前文所示，此处省略以节省篇幅

// 题意概括: 求解方程 B(x)^2 - B(x) + A(x) = 0 (mod x^deg)
// B[0] 由 A[0] 决定: B[0]^2 - B[0] + A[0] = 0.
poly solveQuadFE(poly pa, int deg, ll b0) {
    if (deg == 0) return {};
    poly pb = {b0}; 
    int kl = 1;

    while (kl < deg) {
        kl <<= 1;
        
        poly ca = pa; 
        redeg(ca, kl); 
        
        poly pb2 = mul(pb, pb);
        redeg(pb2, kl);
        
        poly num = sub(pb2, pb); // pb^2 - pb
        num = add(num, ca);       // pb^2 - pb + ca
        redeg(num, kl);

        poly den = smul(pb, 2);  // 2*pb
        if (den.empty() && kl > 0) den.resize(kl,0); else if(den.empty()) den.push_back(0); // 处理pb为空或很短
        if (!den.empty()) den[0] = (den[0] - 1 + MOD) % MOD; // 2*pb - 1
        else { // pb was empty, so 2*pb-1 is -1
            den.assign(kl,0); den[0] = MOD-1;
        }
        redeg(den, kl); 

        poly inv_den = pInv(den, kl); 

        poly term = mul(num, inv_den);
        redeg(term, kl);

        pb = sub(pb, term); 
        redeg(pb, kl);
    }
    redeg(pb, deg);
    return pb;
}

// 主函数示例调用
/*
int main() {
    ios::sync_with_stdio(0); cin.tie(0); cout.tie(0);
    // 预计算itbl等可能需要的全局资源
    pre_itbl(100005); // 假设最大deg是10^5

    int n_deg; // 需要的项数
    cin >> n_deg;

    poly poly_x = {0, 1}; // A(x) = x
    // 对于 T(x)^2 - T(x) + x = 0, T(0)=0 (空树或大小为0的结构计数为0)
    // 或者 T(0)=1 (C(x)=1+xC(x)^2, C(0)=1)
    // 题目是 T(x)^2-T(x)+x=0, T(0)=0.
    // T(x) = x / (1-T(x)) => T - T^2 = x => T^2 - T + x = 0
    // t_0 = 0.
    // t_1 - t_1^2 = x[coeff_x] ...
    // (t_0+t_1x+...)^2 - (t_0+t_1x+...) + x = 0
    // t_0^2 - t_0 = 0 => t_0 = 0 or 1.
    // If t_0=0: (t_1x+t_2x^2)^2 - (t_1x+t_2x^2) + x = 0
    // (t_1^2 x^2 + ...) - t_1x - t_2x^2 + x = 0
    // Coeff of x: -t_1 + 1 = 0 => t_1 = 1.
    // So B[0] = 0.
    poly res_t = solveQuadFE(poly_x, n_deg, 0);

    for (int i = 0; i < n_deg; ++i) {
        cout << (i < (int)res_t.size() ? res_t[i] : 0) << (i == n_deg - 1 ? "" : " ");
    }
    cout << endl;

    return 0;
}
*/
// 需要将所有依赖函数 (qpow, inv, ntt, pInv etc.) 粘贴到此才能完整编译

这个例子展示了如何将一个组合问题转化为函数方程，再利用牛顿迭代求解。

例2：求解

题意概括: 给定一个常系数多项式 和一个多项式 。求解函数方程 。通常 可以从 解出。

解法:
令 。
。
迭代公式：

其中 表示将多项式 代入到 中，即 。
同理。
计算 :
如果 比较小，可以一项项计算 并加权求和。
, etc.
。
这需要 次多项式乘法。如果 是 项，则计算 的复杂度是 。
的计算类似。
总复杂度为 。

Base Case: .
.
这是一个关于 的代数方程。在竞赛中， 通常很小 (如0,1,2)，或者方程有特殊性质使得 容易确定 (例如 是一个解)。
分母 的常数项为 。这必须不为0。

// 辅助函数：计算 P(polyB) = sum c_i * polyB^i
// polyP_coeffs 是 P(y) = sum c_i y^i 的系数 [c0, c1, ..., cD]
poly evalPolyOnPoly(const vector<ll>& polyP_coeffs, const poly& polyB, int deg) {
    poly res(deg, 0);
    poly cur_B_pow = {1}; // B^0 = 1
    redeg(cur_B_pow, deg);

    for (int i = 0; i < (int)polyP_coeffs.size(); ++i) {
        if (i > 0) {
            cur_B_pow = mul(cur_B_pow, polyB); // B^i = B^{i-1} * B
            redeg(cur_B_pow, deg);
        }
        poly term = smul(cur_B_pow, polyP_coeffs[i]); // c_i * B^i
        redeg(term, deg);
        res = add(res, term);
        redeg(res, deg); // 可能不需要，add内部会处理
    }
    return res;
}

// 题意概括: 求解 B(x) = P(B(x)) + A(x) (mod x^deg)
// P(y) = sum pcoeff[i] * y^i
// b0 是 B(0) 的值，需预先解出 b0 = P(b0) + A(0)
poly solveGeneralFE(const vector<ll>& pcoeff, poly pa, int deg, ll b0) {
    if (deg == 0) return {};
    poly pb = {b0};
    int kl = 1;

    // 计算 P'(y) 的系数
    vector<ll> pprime_coeff;
    if (pcoeff.size() > 1) {
        for (size_t i = 1; i < pcoeff.size(); ++i) {
            pprime_coeff.push_back(i * pcoeff[i] % MOD);
        }
    } else { // P(y)是常数，P'(y)=0
        // pprime_coeff remains empty, or {0}
    }


    while (kl < deg) {
        kl <<= 1;
        
        poly ca = pa; redeg(ca, kl); // A(x) mod x^kl
        redeg(pb, kl/2); // pb 是 B_{k-1}, 长度 kl/2
                         // 但在evalPolyOnPoly中作为参数，其项数应为 kl/2
                         // 而evalPolyOnPoly的结果，以及所有中间项应 mod x^kl

        poly p_of_b = evalPolyOnPoly(pcoeff, pb, kl);     // P(B_{k-1}) mod x^kl
        poly pprime_of_b = evalPolyOnPoly(pprime_coeff, pb, kl); // P'(B_{k-1}) mod x^kl

        // Numerator: pb - P(B_{k-1}) - A(x)
        poly num = sub(pb, p_of_b);
        num = sub(num, ca);
        redeg(num, kl);

        // Denominator: 1 - P'(B_{k-1})
        poly den = smul(pprime_of_b, MOD - 1); // -P'(B_{k-1})
        if (den.empty() && kl > 0) den.resize(kl,0); else if(den.empty()) den.push_back(0);
        if (!den.empty()) den[0] = (den[0] + 1 + MOD) % MOD; // 1 - P'(B_{k-1})
        else { den.assign(kl,0); den[0] = 1;}
        redeg(den, kl);
        
        poly inv_den = pInv(den, kl);

        poly term = mul(num, inv_den);
        redeg(term, kl);

        pb = sub(pb, term); // B_k = B_{k-1} - term
        redeg(pb, kl); // 保证 pb 的长度为 kl
    }
    redeg(pb, deg);
    return pb;
}
// 复杂度: T(N) = T(N/2) + O(D * N log N) => O(D * N log N)
// 其中 D 是多项式 P 的度数。
// 空间复杂度: O(N)

注意: 在 evalPolyOnPoly 中，polyB (即迭代中的 ) 的长度是 (上一轮迭代的结果)。但是 evalPolyOnPoly 的结果以及 cur_B_pow 的各项操作都应该在 下进行，所以 redeg 的参数是 kl。这意味着 polyB 虽然只有 项，但在乘法时，比如 ，这个 是 项的，结果 是 项，然后截断到 。

十一、写在最后

多项式牛顿迭代法是一种强大而灵活的工具，它将数值分析中的经典方法巧妙地应用于多项式环，为解决一系列复杂的多项式函数方程提供了统一而高效的框架。通过本文的介绍，我们看到了它在多项式求逆、开方、对数、指数等基本运算中的核心作用，以及它在处理生成函数相关的组合计数问题时的威力。