2021-07-08 23:29 已编辑西北工业大学全栈开发

关注

【博客】浅谈模意义下形式多项式的几种操作

前前言

这篇文章其实是我前年六月份写的，里头仅包含形式多项式的几种基本的操作与部分实现，实际上算法竞赛涉及形式多项式的题目有些已经离谱到十分夸张的程度。

本文充其量只能为接触这方面不久的同学们指一条路，本来我也不准备再搬运到牛客，但无意间看见了牛客另一篇生成函数入门的文章，一方面本文恰巧可以作为那篇文章的后续，另一方面也算为了方便同学们的学习，我还是将这篇文章稍作修订后搬运到了这里。

前言

近几年信息学竞赛中出现了一类利用 生成函数 进行计数的题目，这类题目往往需要一顿推导之后得出一个或几个容易计算的生成函数与代表答案的生成函数间的关系，然后利用这些关系，通过一系列形式多项式的各种操作完成计算。

由于这类题目往往需要将结果对一个大质数取模，因此我们可以用模意义下形式多项式的一些操作来进行这些生成函数之间的运算。

本文将从简单的离散傅里叶变换开始，逐步带你了解这一类基于离散卷积的形式多项式操作。

你需要知道的几点

多项式
- 由数或字母或两者的乘积组成的代数式叫单项式（单个数或字母也为单项式）
  - 单项式中的常数因子为单项式的系数
  - 单项式中所有字母的指数之和为单项式的次数
- 由若干个单项式相加组成的代数式叫多项式
  - 多项式各项中次数最大的单项式的次数为多项式的次数
形式幂级数
- 定义域为 $\mathbb Z^+$ 或其有限子集的函数叫数列（也就是一些有序的数）
- 将数列的每一项求和得到的函数叫级数，一般所讲的级数大多指无穷级数，即数列有无穷项
- 形如的级数叫幂级数（为常数）
  - 其中 $a_i$ 为幂级数的系数（可以类比多项式的各项系数）
  - 若存在一个非负实数 $r$ 使得幂级数在 $| x - x_0 | < r$ 时趋近于一个确定的值（收敛），在 $| x - x_0 | > r$ 时不趋于任何值（发散），则称 $r$ 为幂级数的收敛半径，否则称该幂级数的收敛半径为 $+\infty$
- 将一个数列的各项作为幂级数的系数，所得到的幂级数不一定是收敛的（收敛半径可能为）。但我们淡化其收敛与发散的性质，认为它是收敛的，则可以得到一个形式上类似幂级数的东西（也就是长得像幂级数），我们将其称为形式幂级数
  - 若一个形式幂级数从某一项开始，后面项的系数全部为 $0$ ，则我们称其为形式多项式
离散卷积
- 令 $f(x)$ 和 $g(x)$ 为定义在 $\mathbb R$ 上的两个可积函数，则称 $\int_{-\infty}^{\infty}{f(\tau)g(x-\tau)d\tau}$ 关于 $x$ 的函数为 $f(x)$ 与 $g(x)$ 的卷积
- 与卷积类似，我们称数列 $c_n=\sum_{i=-\infty}^{\infty}{a_ib_{n-i}}$ 为数列 $a_n$ 与数列 $b_n$ 的离散卷积
- 形式幂级数与幂级数一样可以进行运算，其加减运算即是将对应位的系数相加减，而其乘法运算结果定义为两个形式幂级数系数的离散卷积构成的形式幂级数
- 形式幂级数存在乘法逆元当且仅当 $a_0\ne0$ ，且若存在，则逆元唯一

后面的形式幂级数运算若无特殊说明，全部在 $\bmod{x^{n+1}}$ 意义下进行，也即为形式多项式

记号与约定

我们记以数列 $\{f_n\}$ 为系数的形式幂级数为 $F(x)$ 。
我们记 $[x^n]F(x)$ 为 $F(x)$ 的 $n$ 次项系数。
我们记 $f_ig_i$ 为 $f_i\times g_i$ ，而 $F(x)G(x)$ 为 $F(x)$ 与 $G(x)$ 的离散卷积。

多项式乘法

快速傅里叶变换

我们知道，两个形式幂级数做乘法即为它们系数的离散卷积。

假设我们已知 $F(x)$ 和 $G(x)$ 的系数，我们要求得 $H(x)=F(x)G(x)$

根据定义：

$h_n=\sum_{i=0}^{n}{f_ig_{n-i}}\Rightarrow H(x)=\sum_{i=0}^{\infty}{\sum_{j=0}^{i}{f_jg_{i-j}x^i}}$

很显然，我们求出 $H(x)$ 一项的时间复杂度是 $\Theta(n)$ 的，求出 $n$ 项的时间复杂度是 $\Theta(n^2)$ 的，这肯定是不能接受的。我们来考虑怎样更快速地求得 $H(x)$ 。

上面的运算是基于用系数来表示多项式的，事实上，由于 $n+1$ 个点可以唯一确定一个 $n$ 次多项式，我们可以改用 $n+1$ 个点来描述我们的多项式。并且由于我们并不关心形式幂级数的 $x$ 的取值，我们可以任意选取实数甚至虚数作为 $x$ 的值。

这样做的好处是什么呢？我们发现，两个点值表达式相乘只需要将对应点的值相乘即可，也就是说，两个点值表达式做乘法运算是 $\Theta(n)$ 的！

如果我们可以快速地将系数表达式转化成点值表达式，并快速地将点值表达式转化回来，就可以快速地完成多项式乘法。那么具体如何来做呢？

我们先介绍一类神奇的数——单位根 $\omega$ 。

单位根

$n$ 次单位根是 $n$ 次幂为 $1$ 的复数。

也即， $n$ 次单位根 $\omega_n$ 满足 $\omega_n^n=1$ 。

由欧拉公式有 $\omega_n^k=\cos(k\times\frac{2\pi}{n})+i\sin(k\times\frac{2\pi}{n})$ 。

在后文中，我们假设 $n$ 为 $2$ 的正整数次幂：

于是单位根有这几个性质

$\omega_n^{k+\frac n2}=-\omega_n^k$
$\omega_{2n}^{2k}=\omega_n^k$
$\omega_n^0=\omega_n^n=1$

我们来简单证明一下这几个性质：

$\omega_n^{k+\frac n2}=-\omega_n^k$

$\begin{aligned}\omega_n^{k+\frac n2}&=\cos((k+\frac n2)\frac{2\pi}{n})+i\sin((k+\frac n2)\frac{2\pi}{n}) \\&=\cos(k\times\frac{2\pi}{n}+\pi)+i\sin(k\times\frac{2\pi}{n}+\pi) \\&=-(\cos(k\times\frac{2\pi}{n})+i\sin(k\times\frac{2\pi}{n})) \\&=-\omega_n^k\end{aligned}$

$\omega_{2n}^{2k}=\omega_n^k$

$\begin{aligned}\omega_{2n}^{2k}&=\cos(2k\times\frac{2\pi}{2n})+i\sin(2k\times\frac{2\pi}{2n}) \\&=\cos(k\times\frac{2\pi}{n})+i\sin(k\times\frac{2\pi}{n}) \\&=\omega_n^k\end{aligned}$

有了这两个性质，我们就能使用膜法了！

快速傅里叶变换

我们考虑将 $\omega_n$ 的 $0\sim n-1$ 次幂作为 $x$ 代入多项式 $F(x)$ 并求出点值。

也即求出原多项式的离散傅里叶变换（DFT）

我们先将 $F(x)$ 的各项按次数奇偶性分类：

$\begin{aligned}F(x)&=\sum_{i=0}^{n-1}{f_ix^i} \\&=\sum_{i=0}^{\frac n2-1}{f_{2i}x^{2i}}+\sum_{i=0}^{\frac n2-1}{f_{2i+1}x^{2i+1}} \\\end{aligned}$

令 $F_1(x)=\sum_{i=0}^{\frac n2-1}{f_{2i}x^{2i}}$ ， $F_2(x)=\sum_{i=0}^{\frac n2-1}{f_{2i+1}x^{2i+1}}$ ，则有 $F(x)=F_1(x^2)+xF_2(x^2)$ 。

将 $\omega_n^k$ 代入 $x$ ：

$F(\omega_n^k)=F_1(\omega_n^{2k})+\omega_n^kF_2(\omega_n^{2k})$

再将 $\omega_n^{k+\frac n2}$ 代入 $x$ ：

$\begin{aligned}F(\omega_n^{k+\frac n2})&=F_1(\omega_n^{2k+n})+\omega_n^{k+\frac n2}F_2(\omega_n^{2k+n}) \\&=F_1(\omega_n^{2k}\times\omega_n^n)+(\omega_n^k\times\omega_{2n}^n)F_2(\omega_n^{2k}\times\omega_n^n) \\&=F_1(\omega_n^{2k})-\omega_n^kF_2(\omega_n^{2k})\end{aligned}$

发现了吗？这是膜法！

这两个式子间只有第二部分的符号有区别，也就是说，我们计算第一个式子的同时可以直接得到第二个式子的值！

将这种奇偶分类求一半的操作递归下去，我们就得到了一个可以在 $\Theta(n\log n)$ 的时间复杂度内将多项式从系数表达式转化为以单位根各次幂为自变量的点值表达式的算法。

快速傅里叶逆变换

接下来我们再来考虑将这样的点值表达式转化回系数表达式的方法。

也即离散傅里叶逆变换（IDFT）

我们令 $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$ 为 $\{f_0,f_1,f_2,\cdots,f_{n-1}\}$ 的离散傅里叶变换

假设存在 $\{c_0,c_1,c_1,\cdots,c_{n-1}\}$ 满足 $c_k=\sum_{i=0}^{n-1}{y_i(\omega_n^{-k})^i}$

也即 $Y(x)$ 在 $\omega_n^{-k}$ 处的点值，则有：

$\begin{aligned}c_k&=\sum_{i=0}^{n-1}{y_i(\omega_n^{-k})^i} \\&=\sum_{i=0}^{n-1}{(\sum_{j=0}^{n-1}{f_j(\omega_n^i)^j})(\omega_n^{-k})^i} \\&=\sum_{i=0}^{n-1}{\sum_{j=0}^{n-1}{f_j(\omega_n^j)^i}(\omega_n^{-k})^i} \\&=\sum_{i=0}^{n-1}{\sum_{j=0}^{n-1}{f_j(\omega_{n}^{j-k})^i}} \\&=\sum_{j=0}^{n-1}{f_j\sum_{i=0}^{n-1}{(\omega_n^{j-k})^i}}\end{aligned}$

令 $S(x)=\sum_{i=0}^{n-1}{x^i}$ ，将 $\omega_n^k(k\ne0)$ 代入：

$\because\begin{cases}\begin{aligned}S(\omega_n^k)&=\sum_{i=0}^{n-1}{(\omega_n^k)^i} \\\omega_n^kS(\omega_n^k)&=\sum_{i=1}^{n}{(\omega_n^k)^i} \\\end{aligned}\end{cases}$ $\Rightarrow\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-1$ $\therefore S(\omega_n^k)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=0$

而当 $k=0$ 时显然 $S(\omega_n^0)=n$ 。

接着来看之前的式子：

$c_k=\sum_{j=0}^{n-1}{f_j\sum_{i=0}^{n-1}{(\omega_n^{j-k})^i}}$

式子里面的 $\Sigma$ 当 $j\ne k$ 时值为 $0$ ，而 $j=k$ 时值为 $n$ 。因此有 $c_k=nf_k\Rightarrow f_k=\frac{c_k}{n}$

于是只要对原多项式的离散傅里叶变换，以 $\omega_n^{-k}$ 为单位根做一次快速傅里叶变换，再将结果除以 $n$ 即可转化回原多项式。

代码实现

根据上面的理论可以很容易地写出快速傅里叶变换（FFT）的递归实现，但是一般使用时由于递归实现的常数过大，我们一般使用的是迭代实现的快速傅里叶变换。

【例题】[Luogu3803] 多项式乘法评测记录

#include <cstdio>
#include <cmath>

const int N = 2100010;
const double pi = acos(-1);

struct comp {
    double re, im;
    inline comp operator + (const comp &x) const {
        return (comp){re + x.re, im + x.im};
    }
    inline comp operator - (const comp &x) const {
        return (comp){re - x.re, im - x.im};
    }
    inline comp operator * (const comp &x) const {
        return (comp){re * x.re - im * x.im, re * x.im + im * x.re};
    }
} a[N], b[N];
inline void swap(comp &a, comp &b) {
    comp tmp = a; a = b, b = tmp;
}

int lmt, l, r[N];
inline void getRev(int n) {
    lmt = 1, l = 0;
    while (lmt <= n) lmt <<= 1, ++ l;
    for (int i = 1; i < lmt; ++ i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}

inline void DFT(comp *a) {
    comp wt, w, x, y;
    for (int i = 0; i < lmt; ++ i)
        if (i < r[i]) swap(a[i], a[r[i]]);
    for (int i = 1; i < lmt; i <<= 1) {
        wt = (comp){cos(pi / i), sin(pi / i)};
        for (int j = 0, step = i << 1; j < lmt; j += step) {
            w = (comp){1, 0};
            for (int k = 0; k < i; ++ k, w = w * wt) {
                x = a[j + k], y = w * a[i + j + k];
                a[j + k] = x + y, a[i + j + k] = x - y;
            }
        }
    }
}

inline void IDFT(comp *a) {
    comp wt, w, x, y;
    for (int i = 0; i < lmt; ++ i)
        if (i < r[i]) swap(a[i], a[r[i]]);
    for (int i = 1; i < lmt; i <<= 1) {
        wt = (comp){cos(pi / i), -sin(pi / i)};
        for (int j = 0, step = i << 1; j < lmt; j += step) {
            w = (comp){1, 0};
            for (int k = 0; k < i; ++ k, w = w * wt) {
                x = a[j + k], y = w * a[i + j + k];
                a[j + k] = x + y, a[i + j + k] = x - y;
            }
        }
    }
    for (int i = 0; i < lmt; ++ i)
        a[i].re /= lmt, a[i].im /= lmt;
}

int n, m;

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; ++ i)
        scanf("%lf", &a[i].re);
    for (int i = 0; i <= m; ++ i)
        scanf("%lf", &b[i].re);
    getRev(n + m + 2);
    DFT(a); DFT(b);
    for (int i = 0; i < lmt; ++ i)
        a[i] = a[i] * b[i];
    IDFT(a);
    for (int i = 0; i <= n + m; ++ i)
        printf("%.0lf ", a[i].re + 0.1);
    return 0;
}

优化

三次变两次

用刚刚所讲的内容做多项式乘法，需要使用三次 FFT：两次 DFT 和一次 IDFT。

事实上有一个优化可以优化到只做一次 DFT 和一次 IDFT。

我们把 $G(x)$ 的系数放到 $F(x)$ 的虚部里，然后对 $F(x)$ 进行 DFT，将 $F(x)$ 平方后进行 IDFT，将虚部的系数除以 $2$ 即是答案。

考虑为什么这样是对的：对于 $a$ 和 $b$ ，有 $(a+bi)^2=a^2-b^2+(2abi)$ ，因此虚部系数即为答案两倍。

【例题】[Luogu3803] 多项式乘法评测记录

#include <cstdio>
#include <cmath>

const int N = 2100010;
const double pi = acos(-1);

struct comp {
    double re, im;
    inline comp operator + (const comp &x) const {
        return (comp){re + x.re, im + x.im};
    }
    inline comp operator - (const comp &x) const {
        return (comp){re - x.re, im - x.im};
    }
    inline comp operator * (const comp &x) const {
        return (comp){re * x.re - im * x.im, re * x.im + im * x.re};
    }
} a[N];
inline void swap(comp &a, comp &b) {
    comp tmp = a; a = b, b = tmp;
}

int lmt, l, r[N];
inline void getRev(int n) {
    lmt = 1, l = 0;
    while (lmt <= n) lmt <<= 1, ++ l;
    for (int i = 1; i < lmt; ++ i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}

inline void DFT(comp *a) {
    comp wt, w, x, y;
    for (int i = 0; i < lmt; ++ i)
        if (i < r[i]) swap(a[i], a[r[i]]);
    for (int i = 1; i < lmt; i <<= 1) {
        wt = (comp){cos(pi / i), sin(pi / i)};
        for (int j = 0, step = i << 1; j < lmt; j += step) {
            w = (comp){1, 0};
            for (int k = 0; k < i; ++ k, w = w * wt) {
                x = a[j + k], y = w * a[i + j + k];
                a[j + k] = x + y, a[i + j + k] = x - y;
            }
        }
    }
}

inline void IDFT(comp *a) {
    comp wt, w, x, y;
    for (int i = 0; i < lmt; ++ i)
        if (i < r[i]) swap(a[i], a[r[i]]);
    for (int i = 1; i < lmt; i <<= 1) {
        wt = (comp){cos(pi / i), -sin(pi / i)};
        for (int j = 0, step = i << 1; j < lmt; j += step) {
            w = (comp){1, 0};
            for (int k = 0; k < i; ++ k, w = w * wt) {
                x = a[j + k], y = w * a[i + j + k];
                a[j + k] = x + y, a[i + j + k] = x - y;
            }
        }
    }
    for (int i = 0; i < lmt; ++ i)
        a[i].re /= lmt, a[i].im /= lmt;
}

int n, m;

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; ++ i)
        scanf("%lf", &a[i].re);
    for (int i = 0; i <= m; ++ i)
        scanf("%lf", &a[i].im);
    getRev(n + m + 2);
    DFT(a);
    for (int i = 0; i < lmt; ++ i)
        a[i] = a[i] * a[i];
    IDFT(a);
    for (int i = 0; i <= n + m; ++ i)
        printf("%.0lf ", a[i].im / 2 + 0.1);
    return 0;
}

MTT

一种拆系数优化形式多项式乘法的技巧，可以自行查阅资料了解

快速数论变换

如果是在模意义下做运算，我们还有一种算法可以胜任，那就是快速数论变换（NTT）

考虑到做快速傅里叶变换时我们将单位根代入求出点值，如果模意义下有与单位根性质类似的数，是否可以代替单位根完成一样的工作呢？答案是肯定的，这一类数是——原根 $g$ 。

原根

我们先给出阶的定义：

若 $a\perp p$ 且 $p>1$ ， $a$ 模 $p$ 的阶为满足 $a^n\equiv1\pmod{p}$ 的最小的 $n$ ，记为 $\delta_p(a)$ 。

若 $p$ 为正整数， $a$ 是整数，且 $\delta_p(a)=\varphi(p)$ ，则 $a$ 为模 $p$ 的一个原根。（ $\varphi(n)$ 为小于等于 $n$ 的数中与 $n$ 互质的数的个数）

原根有一些很有意思的性质：

若模 $p$ 存在原根，则其一定有 $\varphi(\varphi(p))$ 个原根。
若 $p$ 为质数， $g$ 为 $p$ 的一个原根，则 $g^i\pmod{p}$ 对于不同的 $i\in[1,p)$ 值都不同

同时原根在模意义下与单位根有相同的性质，因此我们可以使用原根代替单位根在模意义下做快速数论变换。

代码实现

类似快速傅里叶变换，只把里面的单位根相应地换成原根，并注意取模即可。

【例题】[Luogu3803] 多项式乘法评测记录

#include <cstdio>

using i64 = long long;

const int MOD = 998244353;
const int N = 2100010;

inline void swap(int &a, int &b) {
    int tmp = a; a = b, b = tmp;
}

inline int pow(int a, int b) {
    int ans = 1;
    while (b) {
        ans = b & 1 ? (i64)ans * a % MOD : ans;
        a = (i64)a * a % MOD; b >>= 1;
    } return ans;
}

int lmt, l, r[N];
inline void getRev(int n) {
    lmt = 1, l = 0;
    while (lmt <= n) lmt <<= 1, ++ l;
    for (int i = 1; i < lmt; ++ i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
}

inline void DFT(int *a) {
    int wt, w, x, y;
    for (int i = 0; i < lmt; ++ i)
        if (i < r[i]) swap(a[i], a[r[i]]);
    for (int i = 1; i < lmt; i <<= 1) {
        wt = pow(3, (MOD - 1) / (i << 1));
        for (int j = 0, step = i << 1; j < lmt; j += step) {
            w = 1;
            for (int k = 0; k < i; ++ k, w = (i64)w * wt % MOD) {
                x = a[j + k], y = (i64)w * a[i + j + k] % MOD;
                a[j + k] = (x + y) % MOD, a[i + j + k] = (x - y + MOD) % MOD;
            }
        }
    }
}

inline void IDFT(int *a) {
    int wt, w, x, y;
    for (int i = 0; i < lmt; ++ i)
        if (i < r[i]) swap(a[i], a[r[i]]);
    for (int i = 1; i < lmt; i <<= 1) {
        wt = pow(332748118, (MOD - 1) / (i << 1));
        for (int j = 0, step = i << 1; j < lmt; j += step) {
            w = 1;
            for (int k = 0; k < i; ++ k, w = (i64)w * wt % MOD) {
                x = a[j + k], y = (i64)w * a[i + j + k] % MOD;
                a[j + k] = (x + y) % MOD, a[i + j + k] = (x - y + MOD) % MOD;
            }
        }
    }
    int bk = pow(lmt, MOD - 2);
    for (int i = 0; i < lmt; ++ i)
        a[i] = (i64)a[i] * bk % MOD;
}

int n, m, a[N], b[N];

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 0; i <= n; ++ i)
        scanf("%d", a + i);
    for (int i = 0; i <= m; ++ i)
        scanf("%d", b + i);
    getRev(n + m + 2);
    DFT(a); DFT(b);
    for (int i = 0; i < lmt; ++ i)
        a[i] = (i64)a[i] * b[i] % MOD;
    IDFT(a);
    for (int i = 0; i <= n + m; ++ i)
        printf("%d ", a[i]);
    return 0;
}

预处理原根

我们来考虑一个优化：在 NTT 的过程中我们多次使用了原根的各次幂，如果我们能将原根的各次幂预处理出来，而不是每次使用都计算一遍，就能有很大的常数优化！

在多项式乘法中我们只做了 $3$ 次 NTT，所以看上去常数优化并不明显。但是在后面的内容中，随着 NTT 次数的增多，预处理原根的优化效果将愈发明显。

【例题】[Luogu3803] 多项式乘法

#include <cstdio>
#include <algorithm>

using std::reverse;
using i64 = long long;
using u64 = unsigned long long;

const int MOD = 998244353;
const int N = 2100010;

inline int pow(int a, int b) {
    int ans = 1;
    while (b) {
        ans = b & 1 ? (i64)ans * a % MOD : ans;
        a = (i64)a * a % MOD; b >>= 1;
    } return ans;
}

int lmt = 1, r[N], w[N];
inline int getLen(int n) {
    return 1 << (32 - __builtin_clz(n));
}

inline void init(int n) {
    int l = 0;
    while (lmt <= n) lmt <<= 1, ++ l;
    for (int i = 1; i < lmt; ++ i)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    int wn = pow(3, (MOD - 1) >> l);
    w[lmt >> 1] = 1;
    for (int i = (lmt >> 1) + 1; i < lmt; ++ i)
        w[i] = (i64)w[i - 1] * wn % MOD;
    for (int i = (lmt >> 1) - 1; i; -- i)
        w[i] = w[i << 1];
    lmt = l;
}

inline void DFT(int *a, int l) {
    static u64 tmp[N];
    int u = lmt - __builtin_ctz(l), t;
    for (int i = 0; i < l; ++ i)
        tmp[r[i] >> u] = a[i];
    for (int i = 1; i < l; i <<= 1)
        for (int j = 0, step = i << 1; j < l; j += step)
            for (int k = 0; k < i; ++ k) {
                t = tmp[i + j + k] * w[i + k] % MOD;
                tmp[i + j + k] = tmp[j + k] + MOD - t;
                tmp[j + k] += t;
            }
    for (int i = 0; i < l; ++ i)
        a[i] = tmp[i] % MOD;
}

inline void IDFT(int *a, int l) {
    reverse(a + 1, a + l); DFT(a, l);
    int bk = MOD - (MOD - 1) / l;
    for (int i = 0; i < l; ++ i)
        a[i] = (i64)a[i] * bk % MOD;
}

int n, m, a[N], b[N], l;

int main() {
    scanf("%d%d", &n, &m);
    init(n + m + 2);
    for (int i = 0; i <= n; ++ i)
        scanf("%d", a + i);
    for (int i = 0; i <= m; ++ i)
        scanf("%d", b + i);
    l = getLen(n + m + 2);
    DFT(a, l); DFT(b, l);
    for (int i = 0; i < l; ++ i)
        a[i] = (i64)a[i] * b[i] % MOD;
    IDFT(a, l);
    for (int i = 0; i <= n + m; ++ i)
        printf("%d ", a[i]);
    return 0;
}

分治 FFT

【例题】[Luogu4721] 分治 FFT

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $g_n=\sum_{i=1}^{n}{g_{n-i}f_i}$ , $g_0=1$ 。

一类存在前后依赖的问题，可以自行查阅资料了解

多项式求导 & 积分

多项式求导

稍微会一点微积分的应该都知道幂法则：

$\frac{d(x^n)}{dx}=nx^{n-1}$

以及和法则：

$\frac{d(f(x)+g(x))}{dx}=\frac{d(f(x))}{dx}+\frac{d(g(x))}{dx}$

因此多项式求导就是将每一项求导后相加，时间复杂度 $\Theta(n)$ 。

代码实现

inline void getDer(int *a, int *b, int deg) {
    for (int i = 0; i + 1 < deg; ++ i)
        b[i] = (i64)a[i + 1] * (i + 1) % MOD;
    b[deg - 1] = 0;
}

多项式积分

与求导同理，有积分公式：

$\int x^ndx=\frac{x^{n+1}}{n+1}+C$

由不定积分的性质有：

$\int(f(x)+g(x))dx=\int f(x)dx+\int g(x)dx$

于是多项式积分也可以每一项积分后相加，时间复杂度 $\Theta(n)$ 。

代码实现

inline void getInt(int *a, int *b, int deg) {
    for (int i = 1; i < deg; ++ i)
        b[i] = (i64)a[i - 1] * inv[i] % MOD;
    b[0] = 0;
}

多项式求逆

【例题】[Luogu4238] 多项式求逆

我们已经知道，当 $f_0\ne0$ 时 $F(x)$ 存在唯一乘法逆元。接下来我们就探讨一下如何求得一个多项式的乘法逆元。

考虑已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $F(x)G(x)\equiv1\pmod{x^n}$

如果 $F(x)$ 只有常数项，答案显然是常数项的乘法逆元，在此基础上，我们考虑倍增地求出整个 $G(x)$ 。

假设我们已经求出了一个 $G'(x)$ 满足：

$F(x)G'(x)\equiv1\pmod{x^{\lceil\frac{n}{2}\rceil}}$

由于 $F(x)G(x)\equiv1\pmod{x^n}$ ，有：

$(G'(x)-G(x))\equiv0\pmod{x^{\lceil\frac n2\rceil}}$

两边同时平方，有：

$(G'(x)-G(x))^2\equiv0\pmod{x^n}$ $G'(x)^2+G(x)^2-2G'(x)G(x)\equiv0\pmod{x^n}$ $F(x)G'(x)^2+G(x)-2G'(x)\equiv0\pmod{x^n}$ $G(x)\equiv2G'(x)-F(x)G'(x)^2\pmod{x^n}$

于是我们得到了从 $G'(x)$ 推到 $G(x)$ 的方式。

时间复杂度 $T(n)=T(\frac n2)+\Theta(n\log n)$ ，由主定理有 $T(n)=\Theta(n\log n)$ 。

代码实现

void getInv(int *a, int *b, int deg) {
    if (deg == 1)
        b[0] = pow(a[0], MOD - 2);
    else {
        static int tmp[N];
        getInv(a, b, (deg + 1) >> 1);
        int len = getLen(deg << 1);
        for (int i = 0; i < deg; ++ i)
            tmp[i] = a[i];
        for (int i = deg; i < len; ++ i)
            tmp[i] = 0;
        DFT(tmp, len); DFT(b, len);
        for (int i = 0; i < len; ++ i)
            b[i] = (2ll - (i64)b[i] * tmp[i] % MOD + MOD) % MOD * b[i] % MOD;
        IDFT(b, len);
        for (int i = deg; i < len; ++ i)
            b[i] = 0;
    }
}

多项式开根

【例题】[Luogu5205] 多项式开根

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)^2\equiv F(x)\pmod{x^n}$ 。

同样如果 $F(x)$ 只有常数项，答案显然是常数项的二次剩余，在此基础上我们仍然考虑倍增。

假设我们已经求出了一个 $G'(x)$ 满足：

$G'(x)^2\equiv F(x)\pmod{x^{\lceil\frac n2\rceil}}$

则有：

$G'(x)^2-G(x)^2\equiv0\pmod(x^{\lceil\frac n2\rceil})$ $G'(x)^4+G(x)^4-2G'(x)^2G(x)^2\equiv0\pmod{x^n}$ $G'(x)^4+G(x)^4+2G'(x)^2G(x)^2\equiv4G'(x)^2G(x)^2\pmod{x^n}$ $G'(x)^2+G(x)^2\equiv2G'(x)G(x)\pmod{x^n}$ $G(x)\equiv\frac{G'(x)^2+G(x)^2}{2G'(x)}\pod{x^n}$

因此有 $G(x)\equiv\frac{G'(x)^2+F(x)}{2G'(x)}\pmod{x^n}$ ，于是可以倍增，时间复杂度同样是 $\Theta(n\log n)$ 。

代码实现

/// 假设常数项为 1
void getSqrt(int *a, int *b, int deg) {
    if (deg == 1) b[0] = 1;
    else {
        static int tmpA[N], tmpB[N];
        getSqrt(a, b, (deg + 1) >> 1);
        getInv(b, tmpB, deg);
        int len = getLen(deg << 1);
        for (int i = 0; i < deg; ++ i)
            tmpA[i] = a[i];
        for (int i = deg; i < len; ++ i)
            tmpA[i] = 0;
        DFT(tmpA, len); DFT(tmpB, len);
        for (int i = 0; i < len; ++ i)
            tmpB[i] = (i64)tmpB[i] * tmpA[i] % MOD;
        IDFT(tmpB, len);
        for (int i = 0; i < deg; ++ i)
            b[i] = (i64)inv[2] * (b[i] + tmpB[i]) % MOD;
        for (int i = 0; i < len; ++ i)
            tmpB[i] = 0;
    }
}

多项式反三角函数（选学）

【例题】[Luogu5265] 多项式反三角函数

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)\equiv\arcsin(F(x))\pmod{x^n}$ 。

如果你对导数表熟悉，你会想起反三角函数的导函数十分地......单纯！

$\frac{d(\arcsin(x))}{dx}=\frac{1}{\sqrt{1-x^2}}$

于是我们考虑先求出 $G(x)$ 的导数再对其积分：

$G(x)\equiv\int{\frac{F'(x)}{\sqrt{1-F(x)^2}}dx}\pmod{x^n}$

$\arccos$ 和 $\arctan$ 也是同理，时间复杂度 $\Theta(n\log n)$ 。

多项式对数函数

【例题】[Luogu4725] 多项式对数函数

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)\equiv\ln(F(x))\pmod{x^n}$ 。

组合意义

什么？你从来没听说过多项式还可以取对数？你不知道它的意义？别着急。

假设我们有一个 $F(x)$ 满足：

$F(x)=\sum_{i=0}^{\infty}{\frac{G(x)^i}{i!}}$

若 $G(x)$ 的第 $i$ 项系数表示集合大小为 $i$ 的方案数，则 $F(x)$ 的第 $i$ 项表示选出若干个集合，使其大小之和恰好为 $i$ 的方案数。

你会发现 $\sum_{i=0}^{\infty}{\frac{G(x)^i}{i!}}$ 恰好为 $e^{G(x)}$ 的麦克劳林级数，于是有：

$F(x)=e^{G(x)}\Longleftrightarrow G(x)=\ln(F(x))$

这就是多项式对数的组合意义。

计算方法

对两边求导，有：

$G'(x)\equiv\frac{F'(x)}{F(x)}\pmod{x^n}$

于是得到：

$G(x)\equiv\int{\frac{F'(x)}{F(x)}dx}\pmod{x^n}$

只要将 $F(x)$ 分别求导和求逆后相乘，将其结果积分即可得到 $G(x)$ ，时间复杂度 $\Theta(n\log n)$ 。

代码实现

/// 假设常数项为 1
inline void getLn(int *a, int *b, int deg) {
    static int tmp[N];
    getInv(a, tmp, deg);
    getDer(a, b, deg);
    int len = getLen(deg << 1);
    DFT(tmp, len); DFT(b, len);
    for (int i = 0; i < len; ++ i)
        tmp[i] = (i64)tmp[i] * b[i] % MOD;
    IDFT(tmp, len);
    getInt(tmp, b, deg);
    for (int i = deg; i < len; ++ i)
        b[i] = 0;
    for (int i = 0; i < len; ++ i)
        tmp[i] = 0;
}

多项式指数函数

【例题】[Luogu4726] 多项式指数函数

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)\equiv e^{F(x)}\pmod{x^n}$ 。

多项式牛顿迭代

想必大家都知道求函数零点除了二分还有一个著名的方法叫牛顿迭代法，也即，求解 $f(x)=0$ ，可以选取一个 $x_0$ 作为根的近似值，然后令 $x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}$ 。

事实上，牛顿迭代法也可以用在多项式身上！

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $F(G(x))\equiv0\pmod{x^n}$ 。

假设我们已经求出了一个 $G_0(x)$ 满足：

$F(G_0(x))\equiv0\pmod{x^{\lceil\frac n2\rceil}}$

对 $F(x)$ 泰勒展开后取前两项，则有：

$F(G(x))\equiv F(G_0(x))+F'(G_0(x))(G(x)-G_0(x))\pmod{x^n}$

可以解出 $G(x)\equiv G_0(x)-\frac{F(G_0(x))}{F'(G_0(x))}\pmod{x^n}$ ，与原来的牛顿迭代公式几乎一致！

事实上之前的多项式开根如果我们使用多项式牛顿迭代可以得出一样的结果！

计算方法

现在回归正题，由 $G(x)\equiv e^{F(x)}\pmod{x^n}$ ：

$\ln(G(x))-F(x)\equiv0\pmod{x^n}$

我们相当于要求上面的方程的零点，对左边求导，可得 $\frac{1}{G(x)}$ ，代入牛顿迭代公式：

$G(x)\equiv G_0(x)-\frac{\ln(G_0(x))-F(x)}{\frac{1}{G_0(x)}}\pmod{x^n}$ $G(x)\equiv G_0(x)(1-\ln(G_0(x))+F(x))\pmod{x^n}$

像前面一样倍增即可，时间复杂度 $\Theta(n\log n)$ 。

代码实现

/// 假设常数项为 0
void getExp(int *a, int *b, int deg) {
    if (deg == 1) b[0] = 1;
    else {
        static int tmp[N];
        getExp(a, b, (deg + 1) >> 1);
        getLn(b, tmp, deg);
        int len = getLen(deg << 1);
        for (int i = 0; i < deg; ++ i)
            tmp[i] = (a[i] - tmp[i] + MOD) % MOD;
        for (int i = deg; i < len; ++ i)
            tmp[i] = 0;
        ++ tmp[0];
        DFT(tmp, len); DFT(b, len);
        for (int i = 0; i < len; ++ i)
            b[i] = (i64)b[i] * tmp[i] % MOD;
        IDFT(b, len);
        for (int i = deg; i < len; ++ i)
            b[i] = tmp[i] = 0;
    }
}

其它

多项式指数函数还有一个基于分治 NTT 的 $\Theta(n\log^2n)$ 做法，有兴趣可以自行查阅资料了解。

多项式双曲函数（选学）

【例题】[Luogu5494] 多项式双曲函数

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)\equiv\sinh(F(x))\pmod{x^n}$ 。

众所周知 $\sinh(x)=\frac{e^x-e^{-x}}{2}$ ，于是我们就可以直接计算：

$G(x)\equiv\frac{e^{F(x)}-e^{-F(x)}}{2}\pmod{x^n}$

$\cosh$ 和 $\operatorname{sech}$ 也是同理，时间复杂度 $\Theta(n\log n)$ 。

多项式三角函数（选学）

【例题】[Luogu5264] 多项式三角函数

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)\equiv\sin(F(x))\pmod{x^n}$ 。

这个看起来好像不是很好做：三角函数的导数仍然是三角函数，而其本身也没有什么容易计算的变形式。

所以，接下来的一切全部都是膜法！

想到三角函数，应该有不少人会想起 Euler 公式：

$e^{i\theta}=\cos(\theta)+i\sin(\theta)$

将 $\theta$ 用 $-\theta$ 替换，可得：

$e^{-i\theta}=\cos(\theta)-i\sin(\theta)$

两式相减，有：

$\begin{aligned}e^{i\theta}-e^{-i\theta}&=2i\sin(\theta) \\\sin(\theta)&=\frac{e^{i\theta}-e^{-i\theta}}{2i}\end{aligned}$

于是有：

$G(x)\equiv\frac{e^{iF(x)}-e^{-iF(x)}}{2i}\pmod{x^n}$

至于 $i$ 的取值，由于 $i^2\equiv-1$ ，解一个二次剩余即可算出来。

同理，也有：

$cos(\theta)=\frac{e^{i\theta}+e^{-i\theta}}{2}$

于是可以计算，时间复杂度 $\Theta(n\log n)$ 。

多项式快速幂

【例题】[Luogu5245] 多项式快速幂

已知多项式 $F(x)$ ，要求出一个多项式 $G(x)$ 满足 $G(x)\equiv F(x)^k\pmod{x^n}$ 。

类似普通快速幂一样做，有一个很容易想的 $\Theta(n\log n\log k)$ 算法。

但是仔细思考一波，发现 $F(x)^k\equiv e^{k\ln(F(x))}\pmod{x^n}$ ，于是问题变得简单。

直接将 $F(x)$ 求对数后系数乘上 $k$ ，再做 $\text{Exp}$ 就是答案，时间复杂度 $\Theta(n\log n)$

多项式高阶前缀和 & 差分

【例题】[Luogu5488] 差分与前缀和

已知多项式 $F(x)$ ，要求出其系数的 $k$ 阶前缀和或差分。

多项式高阶前缀和

考略一阶前缀和的系数：

$[x^i]G_1(x)=\sum_{j=1}^{i}{f_j}$

你会发现这就是 $F(x)$ 和一个系数为全 $1$ 的多项式的离散卷积，即相当于 $G_1(x)=\frac{F(x)}{1-x}$ 。

因此：

$G_k(x)=\frac{F(x)}{(1-x)^k}$

牛顿广义二项式定理告诉我们：

$\begin{aligned}\frac{1}{(1-x)^k}&=\sum_{i=0}^{\infty}{\binom{-k}{i}(-x)^i} \\&=\sum_{i=0}^{n-1}{\binom{i+k-1}{k-1}x^i}\end{aligned}$

于是可以直接计算出这个式子的值，再与 $F(x)$ 进行离散卷积，时间复杂度 $\Theta(n\log n)$ 。

多项式高阶差分

与前缀和相反，一阶差分相当于与 $1-x$ 进行卷积，因此：

$G_k(x)=F(x)(1-x)^k$

使用牛顿二项式定理，有：

$\begin{aligned}(1-x)^k&=\sum_{i=0}^{k}{\binom{k}{i}(-x)^i} \\&=\sum_{i=0}^{k}{(-1)^i\binom{k}{i}x^i}\end{aligned}$

发现这个式子只有 $k+1$ 项，并且可以直接计算，时间复杂度 $\Theta(n\log n)$ 。

多项式除法 & 取模

【例题】[Luogu4512] 多项式除法

已知 $n$ 次多项式 $F(x)$ 和 $m$ 次多项式 $G(x)$ ，要求出多项式 $Q(x),R(x)$ 满足：

$F(x)=Q(x)G(x)+R(x)$

且 $Q(x)$ 为 $n-m$ 次多项式， $R(x)$ 次数小于 $m$ ，保证 $m<n$ 。

我们令 $A_R(x)=x^nA(x^{-1})$ ，可以发现 $A_R(x)$ 即为翻转系数后的 $A(x)$ 。

来考虑题目的式子：

$\begin{aligned}F(x)&=Q(x)G(x)+R(x) \\F(x^{-1})&=Q(x^{-1})G(x^{-1})+R(x^{-1}) \\x^nF(x^{-1})&=x^nQ(x^{-1})G(x^{-1})+x^nR(x^{-1}) \\F_R(x)&=x^{n-m}Q(x^{-1})\cdot x^mG(x^{-1})+x^{n-m+1}\cdot x^{m-1}R(x^{-1}) \\F_R(x)&=Q_R(x)G_R(x)+x^{n-m+1}R_R(x) \\F_R(x)&\equiv Q_R(x)G_R(x)\pmod{x^{n-m+1}} \\Q_R(x)&\equiv F_R(x)G_R(x)^{-1}\pmod{x^{n-m+1}}\end{aligned}$

于是将 $F(x)$ 和 $G(x)$ 系数翻转，让 $F_R(x)$ 与 $G_R(x)$ 的逆元做多项式乘法，结果在 $\bmod{x^{n-m+1}}$ 意义下翻转系数即为 $Q(x)$ 。又由：