多项式

· · 算法·理论

开头先扔板子：多项式板子们

定义

多项式（polynomial）是形如 P(x) = \sum \limits_{i = 0}^{n} a_i x ^ i 的代数表达式。其中 x 是一个不定元。

## 多项式的基本运算 - 多项式的加减法 $$A(x) = \sum_{i = 0}^{n} a_i x ^ i, B(x) = \sum_{i = 0}^{m} b_i x ^ i$$ $$A(x) \pm B(x) = \sum_{i = 0}^{\max(n, m)} (a_i \pm b_i) x ^ i$$ - 多项式的乘法 $$A(x) \times B(x) = \sum_{k = 0}^{nm} \sum_{i + j = k} a_i b_j x ^ k$$ - 多项式除法 这里讨论多项式的带余除法。 可以证明，一定存在唯一的 $q(x), r(x)$ 满足 $A(x) = q(x) B(x) + r(x)$，且 $\partial(r(x)) < \partial(B(x))$。 $q(x)$ 称为 $B(x)$ 除 $A(x)$ 的商式，$r(x)$ 称为 $B(x)$ 除 $A(x)$ 的余式。记作： $$A(x) \equiv r(x) (\bmod \ B(x))$$ 特别的，若 $r(x) = 0$，则称 $B(x)$ 整除 $A(x)$，$A(x)$ 是 $B(x)$ 的一个倍式，$B(x)$ 是 $A(x)$ 的一个因式。记作 $B(x) | A(x)$。 ## 有关多项式的引理 - 对于 $n + 1$ 个点可以唯一确定一个 $n$ 次多项式。 - 如果 $f(x), g(x)$ 都是不超过 $n$ 次的多项式，且它对于 $n +1$ 个不同的数 $\alpha_1, \alpha_2 \cdots \alpha_n$ 有相同的值，即 $\forall i \in [1, n + 1], i \in \mathbb{Z}, f(\alpha_i) = g(\alpha_i)$。则 $f(x) = g(x)$。 ## 多项式的点值表示 如果选取 $n + 1$ 个不同的数 $x_0, x_1, \cdots x_n$ 对多项式进行求值，得到 $A(x_0), A(x_1) \cdots A(x_n)$，那么就称 $\{(x_i, A(x_i)) \ | \ 0 \le i \le n, i \in \mathbb{Z} \}$ 为 $A(x)$ 的点值表示。 ## 快速傅里叶变换（FFT） 快速傅里叶变换是借助单位根进行求值和插值，从而快速进行多项式乘法的算法。 ### 单位根 将复平面上的单位圆均分成 $n$ 份，从 $x$ 轴数，第 $i$ 条线与单位圆的交点称为 $i$ 次单位根，记作 $\omega_{n}^{i}$。 根据定义，可以得到：$\omega_{n}^{1} = i \sin \alpha +\cos \alpha, \alpha = \frac{2 \pi}{n}$。 > 引理：欧拉公式 > $\forall \theta \in R$，$e ^ {i\theta} = i\sin \theta + \cos \theta$。 根据**欧拉公式**，可以得到：$\omega_{n}^{1} = e ^ {\frac{2 \pi i}{n}}$。 由此那么可以得到下面的性质： - $\omega_{n}^{i} \times \omega_{n}^{j} = \omega_{n}^{i + j}

引理：棣莫弗公式

\forall \theta \in R$，$(\cos \theta + i \sin \theta) ^ n = \cos (n\theta) + i \sin (n\theta)

考虑两个复平面上的向量

z_1 = \rho_1(cos \theta_1 + i \sin \theta_1), z_2 = \rho_2(cos \theta_2 + i \sin \theta_2)

计算两个向量相乘的结果，可得：

z_1 z_2 = \rho_1 \rho_2 (\cos (\theta_1 + \theta_2) + i \sin (\theta_1 + \theta_2))

当 \theta_1 = \theta_2, \rho_1 = \rho_2 时便是大名鼎鼎的棣莫弗公式。

因此可以得到，单位圆上的两个复数相乘，模长不变（也就是还在单位圆上），幅角相加。然后可以由此理解上面的三个单位根性质。

离散傅里叶变换（DFT）

离散傅里叶变换，是将 \omega_{n}^{k}, 0 \le k < n 代入 f(x) 和 g(x) 中求值的过程。

对于朴素的方法，每次代入一个单位根，然后用 O(n) 的复杂度计算函数值。时间复杂度 O(n ^ 2)。

离散傅里叶变换利用了单位根的性质巧妙优化了这个过程。离散傅里叶变换过程如下：

首先将 f(x) 根据次数的奇偶性拆成两部分，分别分为：

\begin{cases} e(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i} x ^ {2i} = a_0 + a_2 x^2 + a_4 x^4 \cdots a_{n - 2} x^{n - 2} \\ o(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i + 1} x ^ {2i + 1} = a_1 x + a_3 x^3 + a_5 x^5 \cdots a_{n - 1} x^{n - 1} \end{cases}

设

\begin{cases} e'(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i} x ^ i = a_0 + a_2 x + a_4 x^2 \cdots a_{n - 2} x^{\frac{n - 2}{2}} \\ o'(x) = \sum \limits_{i = 0}^{\frac{n - 2}{2}} a_{2i + 1} x ^ {i} = a_1 + a_3 x^1 + a_5 x^2 \cdots a_{n - 1} x^{\frac{n - 2}{2}} \end{cases}

则 f(x) = e'(x ^ 2) + x o'(x ^ 2)。

将 \omega_{n}^{k} 代入得到：

\begin{cases} f(\omega_{n}^{k}) = e'((\omega_{n}^{k}) ^ 2) + \omega_{n}^{k} o'((\omega_{n}^{k}) ^ 2) = e'(\omega_{n}^{2k}) + \omega_{n}^{k} o'(\omega_{n}^{2k}) \\\\ f(\omega_{n}^{k + \frac{n}{2}}) = e'((\omega_{n}^{k+ \frac{n}{2}}) ^ 2) + \omega_{n}^{k+ \frac{n}{2}} o'((\omega_{n}^{k+ \frac{n}{2}}) ^ 2) = e'(\omega_{n}^{2k}) - \omega_{n}^{k} o'(\omega_{n}^{2k}) \end{cases}

然后你发现，f(\omega_{n}^{k}) 和 f(\omega_{n}^{k + \frac{n}{2}}) 仅仅差了一个符号！！！

所以只需要求出 e'(x) 和 o'(x) 对 \omega^{k}_{n}（0 \le k \le \frac{n}{2}）上的取值，就可以推出 f(x) 在两倍点数上的取值。

每次问题规模缩小一半，因此时间复杂度 O(n \log n)。

离散傅里叶逆变换（IDFT）

假设对于两个多项式都得到了 t = n + m - 1 个点值，设为 \{(x_i, A(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}, \{(x_i, B(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}。

那么可以知道，多项式 C(x) = A(x) \times B(x) 的点值表示一定为：

\{(x_i, A(x_i) B(x_i)) \ | \ 0 \le i < t, i \in \mathbb{Z} \}

现在，只需要将这 t 个点插值回去，就可以得到 A(x)B(x) 了。

先设这 t 个点值分别是：\{(x_i, v_i) \ | \ 0 \le i < t, i \in \mathbb{Z} \}，设最后的多项式为 C(x) = \sum \limits_{i = 0}^{n + m - 2} c_i x^i，这里直接给出结论：

c_k = \dfrac{1}{n} \sum_{i = 0}^{t - 1} v_i \omega_{t}^{-ki}

由此可见，IDFT 和 DFT 仅仅有一个负号的区别。只要将所有的单位根从 \omega_{n}^{k} 变成 \omega_{n}^{-k} 即可。

void FFT(cp a[], int n, int op) {
    if (n == 1) return;
    cp a1[n], a2[n];
    rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
    FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
    cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
    cp Wk = {1, 0};
    rop(i, 0, n >> 1) {
        a[i] = a1[i] + Wk * a2[i];
        a[i + (n >> 1)] = a1[i] - Wk * a2[i];
        Wk = Wk * Wn;
    }
}
void FFT(cp a[], cp b[], int n, int m) {
    m = n + m; n = 1;
    while (n <= m) n <<= 1;
    FFT(a, n, 1); FFT(b, n, 1);
    rop(i, 0, n) a[i] = a[i] * b[i];
    FFT(a, n, -1);
    rep(i, 0, m) a[i].x = a[i].x / n;
}

FFT 优化

void FFT(cp a[], int n, int op) {
    if (n == 1) return;
    cp a1[n], a2[n];
    rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
    FFT(a1, n >> 1, op), FFT(a2, n >> 1, op);
    cp Wn = {cos(2 * pi / n), op * sin(2 * pi / n)};
    cp Wk = {1, 0};
    rop(i, 0, n >> 1) {
        a[i] = a1[i] + a2[i] * Wk;
        a[i + (n >> 1)] = a1[i] - a2[i] * Wk;
        Wk = Wk * Wn;
    }
}
int main() {
    read(n, m);
    rep(i, 0, n) scanf("%lf", &a[i].x);
    rep(i, 0, m) scanf("%lf", &a[i].y);
    m = n + m; n = 1;
    while (n <= m) n <<= 1;
    FFT(a, n, 1);
    rop(i, 0, n) a[i] = a[i] * a[i];
    FFT(a, n, -1);
    rep(i, 0, m) printf("%d ", (int)(a[i].y / 2 / n + 0.5));
}

FFT 例题

  1. A * B problem(plus)

可以设 x = 10，把 a 写成 A(x) = \sum \limits_{i = 0}^{n} a_i x^i 的形式（n = \log_{10} a）。同理可以把 b 转化为多项式 B(x)。

这样求两个数相乘就是求 A(x) \times B(x) 啊。

所以直接 O(n \log n) 做完了。

  1. P3338 [ZJOI2014] 力

给出 n 个数 q_1,q_2, \dots q_n，定义

F_j~=~\sum_{i = 1}^{j - 1} \frac{q_i \times q_j}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i \times q_j}{(i - j)^2} E_i~=~\frac{F_i}{q_i}

对 1 \leq i \leq n，求 E_i 的值。

首先发现这个除以 q_i 就是没用。所以可以化简成：

E_j = \sum_{i = 1}^{j - 1} \frac{q_i}{(i - j)^2}~-~\sum_{i = j + 1}^{n} \frac{q_i}{(i - j)^2}

先看前面这个式子。答案就是：

(j - 1) ^ 2 q_1 + (j - 2) ^ 2 q_2 + (j - 3) ^ 2 q_3 \cdots

设 f(x) = \sum q_i x ^ i, g(x) = \dfrac{1}{i ^ 2} x ^ i。可以发现，E_j' = (f \times g)[j]

再看后面这一块的式子。我们把 f(x) 的系数翻转，变成 f'(x) = \sum q_{n - j + 1} x ^ j。那么可以发现 E_{j}'' = (f' \times g)[n - j + 1]。

跑两次 FFT 就完事了。

  1. P3723 [AH2017/HNOI2017] 礼物

首先发现加减相对于两个手环是对称的。因此可以把对一个手环的减法转化成对另一个手环的加法。这样可以假设全是在第一个手环上执行的加减操作。

第一个手环执行了加 c 的操作，且旋转过之后的序列为 [x_1, x_2 \cdots x_n]，第二个手环为 [y_1, y_2 \cdots y_n]。计算差异值并化简，可以得到差异值是：

\sum x ^ 2 + \sum y ^ 2 + c ^ 2 n + 2c(\sum x - \sum y) - 2 \sum xy

可以发现，这个序列只有最后一项是不定的。

因此将 y 序列翻转后再复制一倍，与 x 卷积，答案就是卷积后序列的 n + 1 \sim 2n 项系数的 \max。

直接暴力枚举 c，加上前面依托就行了。

\text{AC code}

快速数论变换（NTT）

快速数论变换就是基于原根的快速傅里叶变换。

首先考虑快速傅里叶变换用到了单位根的什么性质。

数论中，原根恰好满足这些性质。

对于一个素数的原根 g，设 g_n = g ^ {\frac{p - 1}{n}}。那么：

g_n ^ n = g ^ {p - 1} \equiv 1(\bmod \ p) g_n ^ {\frac{n}{2}} = g ^ {\frac{p - 1}{2}} \equiv - 1(\bmod ~ p) g ^ \alpha + g ^ \beta = g ^ {\alpha + \beta} g_{\alpha n}^{\alpha k} = g_{n}^{k}

我们发现它满足 \omega_{n}^{k} 的全部性质！

因此，只需要用 g_{n}^k = g_{}^{\frac{p - 1}{n} k} 带替全部的 \omega_{n}^{k} 就可以了。

tips：对于一个质数，只有当它可以表示成 p = 2 ^ \alpha + 1，且需要满足多项式的项数 < \alpha 时才能使用 NTT。p 后面有个加一，是因为 g_n 指数的分子上出现了 -1；p - 1 需要时二的整数次幂，是因为每次都要除以 2。

bonus：常用质数及原根

p = 998244353 = 119 \times 2 ^ {23} + 1, g = 3 p = 1004535809 = 479 \times 2 ^ {21} + 1, g = 3 
void NTT(int *a, int n, int op) {
    if (n == 1) return;
    int a1[n], a2[n];
    rop(i, 0, n >> 1) a1[i] = a[i << 1], a2[i] = a[(i << 1) + 1];
    NTT(a1, n >> 1, op), NTT(a2, n >> 1, op);
    int gn = qpow((op == 1 ? g : invg), (mod - 1) / n), gk = 1;
    rop(i, 0, n >> 1) {
        a[i] = (a1[i] + 1ll * gk * a2[i]) % mod;
        a[i + (n >> 1)] = (a1[i] - 1ll * gk * a2[i] % mod + mod) % mod;
        gk = 1ll * gk * gn % mod;
    }
}

NTT 优化：蝴蝶变换

盗大佬一张图

这是进行 NTT 的过程中数组下标的变化。

这样看似乎毫无规律。但是把他们写成二进制：

变换前：000 ~ 001 ~ 010 ~ 011 ~ 100 ~ 101 ~ 110 ~ 111

变换后：000 ~ 100 ~ 010 ~ 110 ~ 001 ~ 101 ~ 011 ~ 111

可以发现，就是对每个下标进行了二进制翻转。

因此可以先预处理出每个下标在递归底层对应的新下标。然后从底层往顶层迭代合并。

预处理每个下标的二进制翻转：

rop(i, 0, n) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit;

优化后的 NTT：

void NTT(int *a, int n, int op) {
    rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
        int gn = qpow((op == 1 ? g : invg), (mod - 1) / (mid << 1));
        for (int i = 0, gk = 1; i < n; i += (mid << 1), gk = 1)
            for (int j = 0; j < mid; j ++ , gk = 1ll * gk * gn % mod) {
                int x = a[i + j], y = 1ll * a[i + j + mid] * gk % mod;
                a[i + j] = Mod(x + y), a[i + j + mid] = Mod(x - y);
            }
    }
}

当然了，FFT 也可以用蝴蝶变换来优化。实践证明，速度变快了至少二分之一。

FFT 的迭代实现

void FFT(cp *a, int n, int op) {
    rop(i, 0, n) if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int mid = 1; (mid << 1) <= n; mid <<= 1) {
        cp Wn = {cos(pi / mid), op * sin(pi / mid)};
        for (int i = 0; i < n; i += (mid << 1)) {
            cp Wk = {1, 0};
            for (int j = 0; j < mid; j ++ , Wk = Wk * Wn) {
                cp x = a[i + j], y = Wk * a[i + j + mid];
                a[i + j] = x + y, a[i + j + mid] = x - y;
            }
        }
    }
}

任意模数多项式乘法（MTT）

计算 f \times g(\bmod ~ p) 的结果（p 不一定能够表示成 2 ^ \alpha + 1 的形式）。

这个东西有两种做法，但是我只学会了三模 NTT。

首先，卷积之后每个系数最多达到 \max \{V\} ^ 2 \times n 的大小。拿模板题举例，这个东西就是 10 ^ {23}。因此只需要找三个模数 p_1, p_2, p_3，满足 p_1p_2p_3 > \max \{V\} ^ 2 \times n，然后用他们分别 NTT，最后再 CRT / exCRT 合并即可。

int CRT(int a, int b, int c, int p) {
    int k = 1ll * Mod(b - a, p2) * inv1 % p2;
    LL x = 1ll * k * p1 + a;
    k = 1ll * Mod((c - x) % p3, p3) * inv2 % p3;
    x = (x + 1ll * k * p1 % p * p2 % p) % p;
    return x;
}
void MTT(int *a, int n, int *b, int m, int p) {
    int B = ((n + m) << 2) + 1;
    rev = new int [B]();
    int *a1 = new int [B](), *b1 = new int [B]();
    int *a2 = new int [B](), *b2 = new int [B]();
    int *a3 = new int [B](), *b3 = new int [B]();
    rop(i, 0, B) 
        a1[i] = a2[i] = a3[i] = a[i], b1[i] = b2[i] = b3[i] = b[i];
    NTT(a1, n, b1, m, p1); NTT(a2, n, b2, m, p2); NTT(a3, n, b3, m, p3);
    inv1 = inv(p1, p2); inv2 = inv(1ll * p1 * p2 % p3, p3);
    rep(i, 0, n + m) a[i] = CRT(a1[i], a2[i], a3[i], p);
}

这里选择的模数为：p_1 = 998244353, p_2 = 469762049, p_3 = 1004535809。他们的原根都为 g = 3。

多项式求逆

求多项式 f(x) 的逆元 f^{-1}(x)。f(x) 的逆元定义为满足 f(x) g(x) = 1(\bmod \ x ^ n) 的多项式 g(x)。

使用倍增法即可求出多项式的逆元。

首先假设已经求出了 f(x) g'(x) \equiv 1(\bmod \ x ^ {\lceil \frac{n}{2} \rceil})。再假设 f(x) \bmod \ x ^ n 意义下逆元为 g(x)。那么有：

g(x) \equiv g'(x) (\bmod \ x ^ {\lceil \frac{n}{2} \rceil}) (g(x) - g'(x)) \equiv 0 (\bmod \ x ^ {\lceil \frac{n}{2} \rceil}) (g(x) - g'(x)) ^ 2 \equiv 0 (\bmod \ x ^ n) g^2(x) - 2 g(x) g'(x) + g'^2(x) \equiv 0 (\bmod \ x ^ {n})

两边同时乘以 f(x)，可以得到：

g(x) - 2 g'(x) + f(x) g'^2(x) \equiv 0 (\bmod \ x ^ n) g(x) \equiv 2g'(x) - f(x) g'^2(x) (\bmod \ x ^ n) g(x) \equiv g'(x) (2 - f(x)g'(x)) (\bmod \ x ^ n)

倍增求即可。

void Inv(int *f, int *g, int n) {
    if (n == 1) {
        g[0] = inv(f[0]); return;
    } Inv(f, g, (n + 1) >> 1);
    int m = 1, bit = 0;
    while (m < (n << 1)) m <<= 1, bit ++ ; bit -- ;
    rop(i, 0, n) tmp[i] = f[i];
    rop(i, n, m) tmp[i] = 0;
    rev = new int [m + 5]();
    rop(i, 1, m) rev[i] = (rev[i >> 1] >> 1) | (i & 1) << bit;
    NTT(tmp, m, 1); NTT(g, m, 1);
    rop(i, 0, m) g[i] = (2ll - 1ll * g[i] * tmp[i] % p + p) % p * g[i] % p;
    NTT(g, m, -1); int invn = inv(m);
    rop(i, 0, m) g[i] = 1ll * g[i] * invn % p;
    rop(i, n, m) g[i] = 0; 
    free(rev); rev = NULL;
}

分治 FFT

给定序列 g_{1\dots n - 1}，求序列 f_{0\dots n - 1}。

其中 f_i=\sum_{j=1}^if_{i-j}g_j，边界为 f_0=1。

答案对 998244353 取模。

其实这是个多项式求逆板子

首先考虑生成函数 F(x) = \sum _{i = 0}^{+ \infty} f_i x ^ i, G(x) = \sum _{i = 0}^{ + \infty}g_i x ^ i。然后可以发现：

F(x) G(x) = \sum \limits_{k = 0}^{+\infty} x ^ k \sum \limits_{i + j = k}^{} f_i g_j = F(x) - f_0 x^0

因此 F(x)(G(x) - 1) = -f_0，也就是：

F(x) \equiv \dfrac{f_0}{1 - G(x)} (\bmod \ x ^ n)

所以直接设 g_0 = 0，然后把 1 - g 求逆就行了。

read(n);
rop(i, 1, n) read(g[i]);
rop(i, 1, n) g[i] = Poly::p - g[i];
(g[0] += 1) %= Poly::p; Poly::Inv(g, n);
rop(i, 0, n) write(' ', g[i]); return 0;

多项式对数函数（Polyln）

给定 f(x)，求多项式 g(x) 使得 g(x) \equiv \ln(f(x)) (\bmod \ x ^ n)

前置知识：简单的求导和积分，以及基本的多项式模板。

首先设 h(x) = \ln(x)，那么

g(x) \equiv h(f(x)) (\bmod \ x ^ n)

然后对同余方程两边同时求导，得到

g'(x) \equiv [h(f(x))]' (\bmod \ x ^ n)

根据复合函数求导公式 f'(g(x)) = f'(g(x))g'(x) 可得，

g'(x) \equiv h'(f(x))f'(x)(\bmod \ x ^ n)

然后又有 \ln '(x) = \dfrac{1}{x}，因此

g'(x) \equiv \dfrac{f'(x)}{f(x)}

计算 f'(x) 和 f^{-1}(x) 作为 a, b。计算 a \times b (\bmod \ 998244353) 得到 g'(x) ，然后求 g'(x) 的积分就好了。

积分公式：\int x ^ a \mathrm{dx} = \dfrac{1}{a + 1} x ^ {a + 1}。

void der(int *f, int n) { rep(i, 1, n) f[i - 1] = 1ll * i * f[i] % p; f[n] = 0; }
void Int(int *f, int n) {dep(i, n, 1) f[i] = 1ll * inv(i) * f[i - 1] % p; f[0] = 0; }
void Ln(int *f, int n) {
    int B = (n << 2) + 5; int *_f = new int[B]();
    rep(i, 0, n) _f[i] = f[i];
    der(_f, n), Inv(f, n);
    NTT(f, n, _f, n); Int(f, n);
    free(_f); _f = NULL;
}

多项式指数函数（PolyExp）

给定多项式 f(x)，求 g(x) 满足 g(x) \equiv e ^ {f(x)} (\bmod \ x ^ n)。

前置知识：牛顿迭代。

牛顿迭代是用来求函数零点的有力工具。

例如，下面这个图是使用牛顿迭代法求 f(x)=x^4+3 x^2+7 x+3 零点的过程。

首先，随便选择一个点 (x_0, f(x_0))，过这个点做 f(x) 的切线。切线方程是 y = f'(x_0)(x - x_0) + f(x_0)。它与 x 轴交于一点 A(-0.56243, 0)。

接下来，再令 x_1 = -0.56243，过点 (x_1, f(x_1)) 再做 f(x) 的切线，与 x 轴交于点 B(-0.60088, 0)。

再令 x_2 = -0.60088，如此迭代下去。可以发现，x_n 会逐渐逼近零点。

刚才说切线方程为 y = f'(x_0)(x - x_0) + f(x_0)。如果令 y = 0，得到的 x 便是切线方程与 x 轴的交点，为

x = x_0 - \dfrac{f(x_0)}{f'(x_0)}

运用于多项式，假设要求使 f(g(x)) \equiv 0(\bmod \ x ^ n) 的多项式 g(x)，并且已经知道 f(g_0(x)) \equiv 0 (\bmod \ x ^ {\frac{n}{2}})。

那么可以说，

g(x) = g_0(x) - \dfrac{f(g_0(x))}{f'(g_0(x))}

接下来解决多项式 Exp。所求为 g(x) 使得 g(x) \equiv e ^ {f(x)} (\bmod \ x ^ n)。两边都取 \ln 得到：

\ln g(x) \equiv f(x) (\bmod \ x ^ n)

移项得：

\ln g(x) - f(x) \equiv 0 (\bmod \ x ^ n)

设 F(g(x)) = ln g(x) - f(x)，那么所求就是 F(x) 的零点。

假设已经有 g_0(x) 使得 F(g_0(x)) \equiv 0 (\bmod \ x ^ {\frac{n}{2}})，根据上面的牛顿迭代，有：

g(x) = g_0(x) - \dfrac{F(g_0(x))}{F'(g_0(x))} = g_0(x) - \dfrac{\ln g_0(x) - f(x)}{\frac{1}{g_0(x)}} = g_0(x)(1 - \ln g_0(x) +f(x))

根据这个柿子倍增求即可。每次需要计算一个 \ln，一个乘法，时间 O(n \log n)。

写的有点丑，超级大肠数。

void Exp(int *f, int *g, int n) {
    if (n == 1) return void(g[0] = 1);
    Exp(f, g, (n + 1) >> 1);
    int B = (n << 2) + 5; int *lnb = new int[B]();
    rop(i, 0, n) lnb[i] = g[i]; Ln(lnb, n);
    tmp = new int[B](); rop(i, 0, n) tmp[i] = f[i];
    rop(i, 0, n) tmp[i] = (1ll * tmp[i] - lnb[i] + p) % p;
    tmp[0] ++ ; NTT(g, n, tmp, n);
    free(tmp); tmp = NULL; free(lnb); lnb = NULL;
}

多项式快速幂（PolyPower）

在模 x ^ n 意义下计算 f^k(x)。

有了前面的知识铺垫，这部分就非常的简单。

根据对数恒等式，有 f(x) = e ^ {\ln f(x)}。

因此 f^k(x) = e ^ {k \ln f(x)}。

把 f(x) 乘以 k，求多项式 \ln，然后再 \exp 回来就行了。

void Power(int *f, int n, int k) {
    Ln(f, n); rop(i, 0, n) f[i] = 1ll * f[i] * k % p; Exp(f, n);
}

多项式开根

在模 x ^ n 意义下计算 f^{\frac{1}{k}}(x)。

这个最简单。直接把 \frac{1}{k} \bmod p，然后多项式快速幂。