【学习笔记】多项式

开个大坑。

FFT 前的前置知识

复数

应该想学FFT的人都会复数吧

复数，即形如 $a+bi$ 的数，其中 $a,b$ 为实数，$i$ 为虚数单位（即 $\sqrt{-1}$）。

在复平面上，一个复数代表着一个点。所谓复平面，即以实数为 $x$ 轴，以虚数为 $y$ 轴所构成的平面直角坐标系。

对于复数加减，直接将实部与虚部相加即可。对于复数乘法，$(a+bi)(c+di)=(ac-bd)+(ad+bc)i$。

在复平面上的复数乘法有这样一个口诀：模相乘，幅角相加。

模即在复平面上点到原点的距离，幅角即相对于 $x$ 轴正半轴所逆时针旋转的角度。

单位根

我们在复平面上做一个半径为 $1$ 的圆，并将他五等分。

这样，我们发现，这些点的模长都为 $1$，并且幅角都是若干倍的关系。

具体的，我们记 $n$ 等分的单位圆上，相对于 $x$ 轴正半轴逆时针旋转的第一个点，记做 ${\omega }_n$，也就是单位根。

根据模相乘，幅角相加，那么第二个点就是 ${\omega }_n^2$，最后一个点就是 ${\omega }_n^{n-1}$。并且，${\omega }_n^n=1$。

单位根的一些性质：

1. ${\omega }_n^k={\omega }_{2n}^{2k}$
   ${\omega }_n^k={\omega }_{\frac{n}{2}}^{\frac{k}{2}}$
   这两个应该比较容易理解。四等分的第二份，相当于八等分的第四份，也相当于两等分的第一份。
2. ${\omega }_{2n}^{n+k}=-{\omega }_{2n}^k$
   因为 ${\omega }_{2n}^n=-1$，所以上式正确。

如何计算单位根？由于是将单位圆进行 $n$ 等分，所以直接使用三角函数计算即可。

即：

$${\omega }_n=\cos \left(\frac{2\pi }{n}\right)+\sin \left(\frac{2\pi }{n}\right)i$$

C++ 中提供了复数类 complex，可以去网上搜索相关用法。不过也可以直接自己写一个，并不难写。

点值表示法

一般，我们所见的多项式的表示方式都是系数表示法，即直接给出 $x^i$ 的系数 $a^i$。而点值表示法，就是给出 $n+1$ 个点，来表示一个多项式。因为给出 $n+1$ 个点后，就可以唯一确定一个 $n$ 次的多项式，参考拉格朗日插值法。

对于多项式乘法来说，系数表示法不好计算，但是点值表示法就很容易计算，直接将相对应的数字相乘即可。所以，我们考虑将系数表示法转换为点值表示法。

离散傅里叶变换

离散傅里叶变换，即 DFT，是一种将系数表示法快速转换为点值表示法的变换。

对于一个 $n$ 次的多项式 $f(x)$：

$$f(x)=a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n$$

我们将它按照奇偶性分开：

$$f(x)=(a_0+a_2{x}^2+a_4{x}^4+\cdots )+(a_{1}x+a_3{x}^3+a_5{x}^5+\cdots )$$

我们分别再设两个多项式 $G(x)=a_0+a_{2}x+a_4{x}^2+\cdots$ 和 $H(x)=a_1+a_{3}x+a_5{x}^2+\cdots$，那么上式就可以表示为：

$$f(x)=G(x^2)+x\cdot H(x^2)$$

我们将 $x={\omega }_n^k$ 代入：

$$\begin{array}{rl}f({\omega }_n^k)& =G(({\omega }_n^k{)}^2)+{\omega }_n^k\cdot H(({\omega }_n^k{)}^2)\\ & =G({\omega }_n^{2k})+{\omega }_n^k\cdot H({\omega }_n^{2k})\\ & =G({\omega }_{\frac{n}{2}}^k)+{\omega }_n^k\cdot H({\omega }_{\frac{n}{2}}^k)\end{array}$$

同理：

$$\begin{array}{rl}f({\omega }_n^{k+\frac{n}{2}})& =G(({\omega }_n^{k+\frac{n}{2}}{)}^2)+{\omega }_n^{k+\frac{n}{2}}\cdot H(({\omega }_n^{k+\frac{n}{2}}{)}^2)\\ & =G((-{\omega }_n^k{)}^2)-{\omega }_n^k\cdot H((-{\omega }_n^k{)}^2)\\ & =G({\omega }_n^{2k})-{\omega }_n^k\cdot H({\omega }_n^{2k})\\ & =G({\omega }_{\frac{n}{2}}^k)-{\omega }_n^k\cdot H({\omega }_{\frac{n}{2}}^k)\end{array}$$

那么我们发现，只要求出 $G({\omega }_{\frac{n}{2}}^k)$ 和 $H({\omega }_{\frac{n}{2}}^k)$ 后，我们就可以同时求出 $f({\omega }_n^{k+\frac{n}{2}})$ 和 $f({\omega }_n^k)$。因此，我们可以直接分治对它进行递归计算。

因为我们进行分治，为了保证左右区间长度相等（不然不好合并），我们可以提前将多项式补到 $2^m-1$ 次，再进行 DFT。

贴一段 OI-Wiki 的代码（我没写过递归版本的）：

#include <cmath>
#include <complex>

typedef std::complex<double> Comp;  // STL complex

const Comp I(0, 1);  // i
const int MAX_N = 1 << 20;

Comp tmp[MAX_N];

void DFT(Comp *f, int n, int rev) {  // rev=1,DFT; rev=-1,IDFT
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  for (int i = 0; i < n; ++i) {  // 偶数放左边，奇数放右边
    if (i & 1)
      f[n / 2 + i / 2] = tmp[i];
    else
      f[i / 2] = tmp[i];
  }
  Comp *g = f, *h = f + n / 2;
  DFT(g, n / 2, rev), DFT(h, n / 2, rev);  // 递归 DFT
  Comp cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI * rev / n));
  // Comp step=exp(I*(2*M_PI/n*rev)); // 两个 step 定义是等价的
  for (int k = 0; k < n / 2; ++k) {
    tmp[k] = g[k] + cur * h[k];
    tmp[k + n / 2] = g[k] - cur * h[k];
    cur *= step;
  }
  for (int i = 0; i < n; ++i) f[i] = tmp[i];
}

时间复杂度是 $O(n\log n)$ 的。

但是这样常数会很大。中间我们不断的对系数进行了交换，我们可以考虑优化这一过程，来观察最后系数位置的变化：

我们将位置化为二进制，来模拟一下这个交换过程：

1. $\{000,001,010,011,100,101,110,111\}$
2. $\{000,010,100,110\}\{001,011,101,111\}$
3. $\{000,100\}\{010,110\}\{001,101\}\{011,111\}$
4. $\{000\}\{100\}\{010\}\{110\}\{001\}\{101\}\{011\}\{111\}$

发现：一个数最后的位置就是一开始的位置的二进制左右翻转。

其实再仔细观察下会发现，第 $i$ 次奇偶性分类其实就是将二进制位 $i$ 为 $0$ 的划分到一块，将 $1$ 划分到一块，原来是高位连续，现在变成了低位连续，于是就是将二进制左右翻转了。

于是，我们可以首先预处理出数组 $r$ 代表二进制位翻转过后的数字。

直接放代码：

for (int i = 0; i < limit; i++)
    r[i] = (r[i >> 1] >> 1) | ((i & 1) * limit >> 1)

自己手模一下就可以，不难理解。

提前交换完毕后，我们发现：合并时，我们将 $k$ 与 $k+\frac{n}{2}$ 合并起来，还是合并到 $k$ 与 $k+\frac{n}{2}$ 这两个位置。所以，我们可以得到 蝴蝶操作：

Complex x = a[k], y = w * a[k + n / 2];
a[k] = x + y, a[k + n / 2] = x - y;

（仅为示例，在实际代码中有所不同）

这样，我们就可以不用任何额外数组，就做到了合并一次答案。最后，我们还可以将递归舍去：直接枚举要合并的区间的长度和区间的位置进行合并。

放代码：

void dft(int limit) {
    for (int i = 0; i < limit; i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 1; mid < limit; mid <<= 1) { // 枚举区间长度
        Complex step(cos(PI / mid), sin(PI / mid)); // 单位根
        for (int l = 0, len = mid << 1; l < limit; l += len) { // 枚举区间左端点
            Complex cur(1, 0);
            for (int k = 0; k < mid; k++, cur = cur * step) { // 枚举 k
                Complex x = a[l + k], y = cur * a[l + k + mid];
                a[l + k] = x + y;
                a[l + k + mid] = x - y;
            }
        }
    }
}

逆离散傅里叶变换

逆离散傅里叶变换，即 IDFT，就是将点值表示法转回系数表示法的方法。

我们从线性代数的角度考虑上面的 DFT 操作，其实就是将系数向量乘上了一个矩阵，也就是：

$$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ y_3\\ ⋮\\ y_{n-1}\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& {\omega }_n^3& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& {\omega }_n^6& \cdots & {\omega }_n^{2(n-1)}\\ 1& {\omega }_n^3& {\omega }_n^6& {\omega }_n^9& \cdots & {\omega }_n^{3(n-1)}\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& {\omega }_n^{3(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ ⋮\\ a_{n-1}\end{array}\right]$$

我们记左面的点值向量为 $E$，中间的矩阵为 $D$，右面的系数矩阵为 $V$，那么 $D\times V=E$。

我们求 $V$，实际上就是要求 $D$ 的逆矩阵。

考虑构造 $A_{ij}=D_{ij}^{-1}$，来计算 $A\times D$：

$$\begin{array}{rl}(A\times D{)}_{ij}& =\sum _{k=0}^{n-1}{A}_{ik}{D}_{kj}\\ & =\sum _{k=0}^{n-1}{\omega }_n^{-ik}{\omega }_n^{kj}\\ & =\sum _{k=0}^{n-1}{\omega }_n^{k(j-i)}\\ & =\{\begin{array}{ll}n& (i=j)\\ \sum _{k=0}^{n-1}({\omega }_n^{j-i}{)}^k=\frac{({\omega }_n^{j-i}{)}^n-1}{{\omega }_n^{j-i}-1}=\frac{({\omega }_n^n{)}^{j-i}-1}{{\omega }_n^{j-i}-1}=0& (i\ne j)\end{array}\end{array}$$

发现这个矩阵就是单位矩阵的 $n$ 倍。即 $\frac{A\times D}{n}=I$，那么 $D^{-1}=\frac{1}{n}A$。

这样，我们将点值向量左乘一个 $\frac{1}{n}A$，就得到了系数表示法。

具体实现上，由于 $A_{ij}=D_{ij}^{-1}={\omega }_n^{-ij}={\omega }_n^{n-ij}$，它与 ${\omega }_n^{ij}$ 在复平面上是关于 $x$ 轴对称的，所以直接对虚部取一个相反数即可。那么，我们将 DFT 函数中引入一个 type 参数，表示是否是逆操作，那么就可以写出下面的代码：

void dft(int limit, int type) {
    for (int i = 0; i < limit; i++) if (i < r[i]) swap(a[i], a[r[i]]);
    for (int mid = 1; mid < limit; mid <<= 1) { // 枚举区间长度
        Complex step(cos(PI / mid), type * sin(PI / mid)); // 单位根
        for (int l = 0, len = mid << 1; l < limit; l += len) { // 枚举区间左端点
            Complex cur(1, 0);
            for (int k = 0; k < mid; k++, cur = cur * step) { // 枚举 k
                Complex x = a[l + k], y = cur * a[l + k + mid];
                a[l + k] = x + y;
                a[l + k + mid] = x - y;
            }
        }
    }
    if (type == -1) // 不要忘记 1/n
        for (int i = 0; i < limit; i++) a[i].r /= limit;
}

附封装版的完整代码：

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 60005;
struct Complex {
    double r, i;
    Complex(double r = 0, double i = 0) : r(r), i(i) {}
    Complex operator+(Complex b) { return { r + b.r, i + b.i }; }
    Complex operator-(Complex b) { return { r - b.r, i - b.i }; }
    Complex operator*(Complex b) { return { r * b.r - i * b.i, r * b.i + i * b.r }; }
};
int r[MAXN];
const double PI = acos(-1.0);
struct Polynomial {
    vector<Complex> a;
    int len;
    Complex& operator[](int b) { return a[b]; }
    Polynomial(int len = 0) : len(len) { a.resize(len + 1); }
    void set(int len) { this->len = len, a.resize(len + 1); }
    void dft(int limit, int type) {
        set(limit);
        for (int i = 0; i < limit; i++) 
            if (i < r[i]) swap(a[i], a[r[i]]);
        for (int mid = 1; mid < limit; mid <<= 1) {
            Complex step(cos(PI / mid), type * sin(PI / mid));
            for (int l = 0; l < limit; l += (mid << 1)) {
                Complex w(1, 0);
                for (int j = 0; j < mid; j++, w = w * step) {
                    Complex x = a[l + j], y = w * a[l + j + mid];
                    a[l + j] = x + y, a[l + j + mid] = x - y;
                }
            }
        }
        if (type == -1) for (int i = 0; i < limit; i++) a[i].r /= limit;
    }
    Polynomial operator*(Polynomial b) {
        Polynomial a = *this, c; 
        int n = len + b.len;
        int limit = 1; 
        while (limit <= n) limit <<= 1; 
        c.set(limit);
        for (int i = 0; i < limit; i++) 
            r[i] = (r[i >> 1] >> 1) | ((i & 1) * limit >> 1);
        a.dft(limit, 1), b.dft(limit, 1);
        for (int i = 0; i < limit; i++) c[i] = a[i] * b[i];
        c.dft(limit, -1);
        c.set(n);
        return c;
    }
};

NTT 快速数论变换

FFT 有一个重要的问题：精度问题。毕竟涉及到浮点数运算，精度误差是不可避免的，而有时候我们需要求取模意义下的数，这时候普通的 FFT 就不适用了。这时候，我们就引出了 NTT。

精度的瓶颈在哪？显然是单位根的运算。我们考虑模意义下什么东西可以代替单位根。

——没错，就是原根。

不明白原根的建议自行百度，上一篇博客写到了原根，但是写的并不详细。

我们来一一看这几条性质：

1. ${\omega }_n^k$ 对于 $k\in [0,n-1]$ 互不相等。
   若模数 $p$ 为一个质数，那么对于原根 $g$，$g^k,k\in [0,p-1]$ 互不相等，那么 $(g^{\frac{p-1}{n}}{)}^k,k\in [0,n-1]$ 也互不相等。
2. ${\omega }_n^k={\omega }_{\frac{n}{2}}^{\frac{k}{2}}$
   $$(g^{\frac{p-1}{n/2}}{)}^{\frac{k}{2}}=g^{\frac{p-1}{n/2}\times \frac{k}{2}}=(g^{\frac{p-1}{n}}{)}^k$$

3. ${\omega }_{2n}^{k+n}=-{\omega }_{2n}^k$
   由费马小定理可得，$g^{p-1}\equiv 1$，那么 $g^{p-1}-1\equiv (g^{\frac{p-1}{2}}-1)(g^{\frac{p-1}{2}}+1)\equiv 0$，即 $g^{\frac{p-1}{2}}=\pm 1$。
   又因为 $g^0\equiv 1$，根据性质 1，$g^{\frac{p-1}{2}}=-1$，那么 $(g^{\frac{p-1}{2n}}{)}^n\equiv g^{\frac{p-1}{2}}\equiv -1$，即 $(g^{\frac{p-1}{2n}}{)}^{n+k}\equiv -(g^{\frac{p-1}{2n}}{)}^k$。

于是，我们可以得出以下结论：原根就是模意义下的单位根！

不过我们发现，$n$ 是 $2^m$ 的形式，那么 $\frac{p-1}{n}$ 要想是整数，那么 $p$ 必须是 $r\times 2^t+1$ 的形式。

实际上，有：

$$998244353=119\times 2^{23}+1,g=3$$

另外一些常用的 NTT 模数：

$$469762049=7\times 2^{26}+1,g=3$$

$$1004535809=479\times 2^{21}+1,g=3$$

考场上如何检验一个数是不是 NTT 模数？

gnome-calculator 有个功能叫质因数分解

如果模数不是 NTT 模数怎么办？

那你就大骂出题人↓↑

任意模数多项式乘法

NTT code:

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 3000005, P = 998244353, G = 3, GI = 332748118;
int qpow(int a, int b) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = 1ll * ans * a % P;
        a = 1ll * a * a % P;
        b >>= 1;
    }
    return ans;
}
int r[MAXN];
struct Polynomial {
    vector<int> a;
    int len;
    int& operator[](int b) { return a[b]; }
    Polynomial(int len = 0) : len(len) { a.resize(len + 1); }
    void set(int b) { len = b, a.resize(b + 1); }
    void ntt(int limit, bool rev) {
        set(limit);
        for (int i = 0; i < limit; i++)
            if (i < r[i]) swap(a[i], a[r[i]]);
        for (int mid = 1; mid < limit; mid <<= 1) {
            int step = qpow(rev ? GI : G, (P - 1) / (mid << 1));
            for (int l = 0; l < limit; l += (mid << 1)) {
                int w = 1;
                for (int j = 0; j < mid; j++, w = 1ll * w * step % P) {
                    int x = a[l + j], y = 1ll * w * a[l + j + mid] % P;
                    a[l + j] = (x + y) % P, a[l + j + mid] = (1ll * P + x - y) % P;
                }
            }
        }
        if (rev) {
            int ninv = qpow(limit, P - 2);
            for (int i = 0; i < limit; i++) 
                a[i] = 1ll * a[i] * ninv % P;
        }
    }
    Polynomial operator*(Polynomial b) {
        Polynomial a = *this, c; int n = a.len + b.len;
        int limit = 1; while (limit <= n) limit <<= 1;
        c.set(limit);
        for (int i = 1; i < limit; i++)
            r[i] = (r[i >> 1] >> 1) | ((i & 1) * limit >> 1);
        a.ntt(limit, 0);
        b.ntt(limit, 0);
        for (int i = 0; i <= limit; i++) c[i] = 1ll * a[i] * b[i] % P;
        c.ntt(limit, 1);
        c.set(n);
        return c;
    }
};

应用

多项式乘法通常可以用于快速计算卷积。

观察多项式乘法的式子：

$$F\times G=\sum _{i=0}^n\sum _{j=0}^m{f}_i{g}_j{x}^{i+j}$$

我们可以转换一下：

$$F\times G=\sum _{i=0}^{n+m}\sum _{j=0}^i{f}_j{g}_{i-j}{x}^i$$

我们发现，$x^i$ 的系数为 $\sum _{j=0}^n{f}_j{g}_{i-j}$，也就是卷积的形式。

于是，我们可以将很多形如卷积形式的式子使用 FFT/NTT 进行优化。

例题：

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\le i\le n$，求 $E_i$​ 的值。

首先化式子：

$$\begin{array}{rl}{E}_i& =\frac{F_i}{q_i}\\ & =\sum _{j=1}^{i-1}\frac{q_j}{(j-i{)}^2}-\sum _{j=i+1}^n\frac{q_j}{(j-i{)}^2}\\ & =\sum _{j=1}^{i-1}\frac{q_j}{(i-j{)}^2}-\sum _{j=1}^{n-i}\frac{q_{i+j}}{j^2}& (j\to j+i)\end{array}$$

左面显然是卷积形式，设 $F(x)=q_x,G(x)=\frac{1}{x^2}$，那么左面的式子就是 $(F\cdot G)(i-1)$。

那么看右面：右面是差一定的形式，而不是和一定。我们可以将 $F(x)$ 的下标翻转过来，即设 $F^{\prime }(x)=F(n-x)$，那么右面就变为了

$$\sum _{j=1}^{n-i}{F}^{\prime }(n-i-j)G(j)$$

发现这就也是卷积的形式了，直接计算即可。

P4199 万径人踪灭

见 上篇文章。

全 家 桶 ！

实际上我就学了几个，为了防止过段时间全忘了，就先记一下。

多项式乘法逆

给定多项式 $F(x)$，求一个多项式 $G(x)$ 使得 $F(x)G(x)\equiv 1(\mathrm{mod}{x}^n)$。

考虑倍增。假如现在已经求出了 $F(x)H(x)\equiv 1(\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })$：

$$\begin{array}{rlr}F(x)H(x)& \equiv 1& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ F(x)G(x)& \equiv 1& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ F(x)(G(x)-H(x))& \equiv 0& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ G(x)-H(x)& \equiv 0& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ (G(x)-H(x){)}^2& \equiv 0& (\mathrm{mod}{x}^n)\\ G(x{)}^2-2G(x)H(x)+H(x{)}^2& \equiv 0& (\mathrm{mod}{x}^n)\\ F(x)G(x{)}^2-2F(x)G(x)H(x)+F(x)H(x{)}^2& \equiv 0& (\mathrm{mod}{x}^n)\\ G(x)-2H(x)+F(x)H(x{)}^2& \equiv 0& (\mathrm{mod}{x}^n)\\ G(x)& \equiv 2H(x)-F(x)H(x{)}^2& (\mathrm{mod}{x}^n)\end{array}$$

$n=1$ 时，直接求逆元即可。

Polynomial inv(int n) {
    Polynomial b;
    b[0] = qpow(a[0], P - 2);
    for (int d = 1; d < (n << 1); d <<= 1) {
        Polynomial a = *this;
        a.set(d - 1);
        int limit = d << 1;
        calcRev(limit);
        a.ntt(limit, 0), b.ntt(limit, 0);
        for (int i = 0; i < limit; i++) 
            b[i] = (2ll - 1ll * a[i] * b[i] % P + P) % P * b[i] % P;
        // 直接对点值进行计算，而不是两次多项式乘法，优化常数
        b.ntt(limit, 1);
        b.set(d - 1);
    }
    b.set(n - 1);
    return b;
}

多项式开根

给定多项式 $F(x)$，求一个多项式 $G(x)$ 使得 $G(x{)}^2\equiv F(x)(\mathrm{mod}{x}^n)$。

同样考虑倍增：假设已经求得 $H(x{)}^2\equiv F(x)(\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })$，那么：

$$\begin{array}{rlr}H(x{)}^2& \equiv F(x)& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ G(x{)}^2& \equiv F(x)& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ H(x{)}^2& \equiv G(x{)}^2& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ H(x)& \equiv G(x)& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ H(x)-G(x)& \equiv 0& (\mathrm{mod}{x}^{\lceil \frac{n}{2}\rceil })\\ (H(x)-G(x){)}^2& \equiv 0& (\mathrm{mod}{x}^n)\\ H(x{)}^2-2H(x)G(x)+G(x{)}^2& \equiv 0& (\mathrm{mod}{x}^n)\\ H(x{)}^2-2H(x)G(x)+F(x)& \equiv 0& (\mathrm{mod}{x}^n)\\ G(x)& \equiv \frac{H(x{)}^2+F(x)}{2H(x)}& (\mathrm{mod}{x}^n)\end{array}$$

求逆元，然后倍增算就可以了。

对于 $n=1$ 的情况，在强化版中，要计算二次剩余 $x^2\equiv a_0$，我不会，咕了。

Polynomial sqrt(int n) {
    static int TWOINV = qpow(2, P - 2);
    Polynomial b;
    b[0] = 1;
    for (int d = 1; d < (n << 1); d <<= 1) {
        Polynomial a = *this, c = b.inv(d);
        a.set(d - 1);
        int limit = d << 1;
        calcRev(limit);
        a.ntt(limit, 0), c.ntt(limit, 0);
        for (int i = 0; i < limit; i++) a[i] = 1ll * a[i] * c[i] % P;
        a.ntt(limit, 1);
        b.set(d - 1);
        for (int i = 0; i < d; i++) b[i] = 1ll * (a[i] + b[i]) * TWOINV % P;
    }
    b.set(n - 1);
    return b;
}

多项式对数函数（多项式 $\ln$）

给定多项式 $F(x)$，求一个多项式 $G(x)$ 使得 $G(x)\equiv \ln F(x)(\mathrm{mod}{x}^n)$。

直接求不好求，我们给他求个导：

$$G^{\prime }(x)\equiv \frac{F^{\prime }(x)}{F(x)}(\mathrm{mod}{x}^n)$$

再积分回去：

$$G(x)\equiv \int \frac{F^{\prime }(x)}{F(x)}(\mathrm{mod}{x}^n)\mathrm{d}x$$

求导与积分：

Polynomial derivative() {
    Polynomial b(len - 1);
    for (int i = 1; i <= len; i++) b[i - 1] = 1ll * a[i] * i % P;
    return b;
}
Polynomial integral() {
    Polynomial b(len + 1);
    for (int i = 0; i <= len; i++) b[i + 1] = 1ll * a[i] * qpow(i + 1, P - 2) % P;
    return b;
}

对数：

Polynomial ln(int n) {
    return (derivative() * inv(n)).integral().set(n - 1);
}

多项式指数函数（多项式 $\exp$）

需要用一个东西：

牛顿迭代：

什么是牛顿迭代？

牛顿迭代是用来求零点的方法。首先有一个近似值 $x_0$，做 $x_0$ 处的切线，交 $x$ 轴于一点，再以这点作为 $x_0$ 继续迭代。

假设函数是 $f(x)$，我们写出切线方程：

$$y=f^{\prime }(x_0)(x-x_0)+f(x_0)$$

令 $y=0$，那么

$$x=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$$

将其运用到多项式上，就是：

$$G(x)=G_0(x)-\frac{F(G_0(x))}{F(G_0(x))}$$

这样每一次迭代精度会翻倍（证明见 OI-Wiki），于是我们可以使用牛顿迭代来进行一些多项式操作。

例如此题：

$$G(x)\equiv e^{F(x)}(\mathrm{mod}{x}^n)$$

我们两边取对数：

$$\begin{array}{rlr}\ln G(x)& \equiv F(x)& (\mathrm{mod}{x}^n)\\ \ln G(x)-F(x)& \equiv 0& (\mathrm{mod}{x}^n)\end{array}$$

于是我们设 $H(G(x))=\ln G(x)-F(x)$，由于 $F(x)$ 在这里是常数，那么求导后得到 $H^{\prime }(G(x))=\frac{1}{G(x)}$，那么牛顿迭代：

$$\begin{array}{rl}G(x)& =G_0(x)-\frac{\ln G_0(x)-F(x)}{\frac{1}{G_0(x)}}\\ & =G_0(x)(1-\ln G_0(x)+F(x))\end{array}$$

于是继续倍增就可以了！

Polynomial exp(int n) {
    Polynomial b;
    b[0] = 1;
    for (int d = 1; d < (n << 1); d <<= 1) {
        Polynomial a = *this, e = b.ln(d);
        a.set(d - 1);
        b = b * (e * (P - 1) + a + 1);
        b.set(d - 1);
    }
    b.set(n - 1);
    return b;
}

多项式除法

给定 $n$ 次多项式 $F(x)$ 和 $m$ 次多项式 $G(x)$，求一个 $n-m$ 次多项式 $Q(x)$ 和一个 $m-1$ 次多项式 $R(x)$，满足 $F(x)=Q(x)G(x)+R(x)$。

我们设 $F_R(x)$ 为将 $F(x)$ 系数翻转之后的式子，也就是：

$$F_R(x)=\sum _{i=0}^n{a}_{n-i}{x}^i$$

我们将 $\frac{1}{x}$ 带入：

$$\begin{array}{rl}F\left(\frac{1}{x}\right)& =\sum _{i=0}^n{a}_i{\left(\frac{1}{x}\right)}^i\\ & =\sum _{i=0}^n{a}_i{x}^{-i}\\ & =\sum _{i=0}^n{a}_{n-i}{x}^{i-n}\\ & =\frac{1}{x^n}\sum _{i=0}^n{a}_{n-i}{x}^i\\ F\left(\frac{1}{x}\right)x^n& =F_R(x)\end{array}$$

那么有：

$$\begin{array}{rl}F\left(\frac{1}{x}\right)& =G\left(\frac{1}{x}\right)Q\left(\frac{1}{x}\right)+R\left(\frac{1}{x}\right)\\ F\left(\frac{1}{x}\right)x^n& =G\left(\frac{1}{x}\right)Q\left(\frac{1}{x}\right)x^n+R\left(\frac{1}{x}\right)x^n\\ F\left(\frac{1}{x}\right)x^n& =G\left(\frac{1}{x}\right)x^{m}Q\left(\frac{1}{x}\right)x^{n-m}+R\left(\frac{1}{x}\right)x^{m-1}\cdot x^{n-m+1}\\ F_R(x)& =G_R(x)Q_R(x)+R_R(x)\cdot x^{n-m+1}\end{array}$$

$$\begin{array}{rl}{F}_R(x)& \equiv G_R(x)Q_R(x)(\mathrm{mod}{x}^{n-m+1})\\ Q_R(x)& \equiv \frac{F_R(x)}{G_R(x)}(\mathrm{mod}{x}^{n-m+1})\end{array}$$

而 $Q(x)$ 正好是 $n-m$ 次的，所以直接计算即可得到 $Q(x)$。

那么 $R(x)=F(x)-G(x)Q(x)$，即可计算出 $R(x)$。

void reverse() { std::reverse(a.begin(), a.end()); }
pair<Polynomial, Polynomial> operator/(Polynomial b) {
    int n = len, m = b.len;
    Polynomial ra = *this, rb = b, rc, d;
    ra.reverse(), rb.reverse();
    rc = ra * rb.inv(n - m + 1);
    rc.set(n - m);
    rc.reverse();
    d = *this - b * rc;
    d.set(m - 1);
    return make_pair(rc, d);
}

附上我目前的模板：

#include <bits/stdc++.h>
using namespace std;
const int MAXN = 1200005, P = 998244353, G = 3, GI = 332748118;
const int INV2 = (P + 1) / 2;
int qpow(int a, int b) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = 1ll * ans * a % P;
        a = 1ll * a * a % P;
        b >>= 1;
    }
    return ans;
}
int r[MAXN];
struct Polynomial {
    vector<int> a;
    int len;
    int& operator[](int b) { return a[b]; }
    Polynomial(int len = 0) : len(len) { a.resize(len + 1); }
    Polynomial(initializer_list<int> l) : len(l.size() - 1), a(l) {}
    Polynomial& set(int b) { len = b, a.resize(b + 1); return *this; }
    static Polynomial resize(Polynomial a, int s) { Polynomial b = a; return b.set(s); }
    void reverse() { std::reverse(a.begin(), a.end()); }
    static void calcRev(int limit) {
        for (int i = 1; i < limit; i++)
            r[i] = (r[i >> 1] >> 1) | ((i & 1) * limit >> 1);
    }
    void ntt(int limit, bool rev) {
        set(limit);
        for (int i = 0; i < limit; i++)
            if (i < r[i]) swap(a[i], a[r[i]]);
        for (int mid = 1; mid < limit; mid <<= 1) {
            int step = qpow(rev ? GI : G, (P - 1) / (mid << 1));
            for (int l = 0; l < limit; l += (mid << 1)) {
                int w = 1;
                for (int j = 0; j < mid; j++, w = 1ll * w * step % P) {
                    int x = a[l + j], y = 1ll * w * a[l + j + mid] % P;
                    a[l + j] = (x + y) % P, a[l + j + mid] = (1ll * P + x - y) % P;
                }
            }
        }
        int invn = qpow(limit, P - 2);
        if (rev) {
            for (int i = 0; i < limit; i++) 
                a[i] = 1ll * a[i] * invn % P;
        }
    }
    void print() { for (int i : a) printf("%d ", i); printf("\n"); }
    Polynomial operator*(Polynomial b) {
        Polynomial a = *this, c;
        int n = a.len + b.len;
        int limit = 1;
        while (limit <= n) limit <<= 1;
        c.set(limit);
        calcRev(limit);
        a.ntt(limit, 0), b.ntt(limit, 0);
        for (int i = 0; i < limit; i++) c[i] = 1ll * a[i] * b[i] % P;
        c.ntt(limit, 1);
        c.set(n);
        return c;
    }
    Polynomial operator*(int b) {
        Polynomial c = *this;
        for (int& i : c.a) i = 1ll * i * b % P;
        return c;
    }
    Polynomial operator+(int b) {
        Polynomial c = *this;
        c[0] = (1ll * c[0] + b + P) % P;
        return c;
    }
    Polynomial operator-(int b) {
        Polynomial c = *this;
        c[0] = (1ll * c[0] - b + P) % P;
        return c;
    }
    Polynomial operator+(Polynomial b) {
        Polynomial c = *this;
        c.set(max(c.len, b.len));
        for (int i = 0; i <= b.len; i++) c[i] = (c[i] + b[i]) % P;
        return c;
    }
    Polynomial operator-(Polynomial b) {
        Polynomial c = *this;
        c.set(max(c.len, b.len));
        for (int i = 0; i <= b.len; i++) c[i] = (c[i] - b[i] + P) % P;
        return c;
    }
    Polynomial shift(int x) {
        Polynomial b(len - x);
        for (int i = max(x, 0); i <= len; i++)
            b[i - x] = a[i];
        return b;
    }
    Polynomial inv(int n) {
        Polynomial b;
        b[0] = qpow(a[0], P - 2);
        for (int d = 1; d < (n << 1); d <<= 1) {
            Polynomial a = *this;
            a.set(d - 1);
            int limit = d << 1;
            calcRev(limit);
            a.ntt(limit, 0), b.ntt(limit, 0);
            for (int i = 0; i < limit; i++) b[i] = (2ll - 1ll * a[i] * b[i] % P + P) % P * b[i] % P;
            b.ntt(limit, 1);
            b.set(d - 1);
        }
        b.set(n - 1);
        return b;
    }
    Polynomial sqrt(int n) {
        static int INV2 = qpow(2, P - 2);
        Polynomial b;
        b[0] = 1;
        for (int d = 1; d < (n << 1); d <<= 1) {
            Polynomial a = *this, c = b.inv(d);
            a.set(d - 1);
            int limit = d << 1;
            calcRev(limit);
            a.ntt(limit, 0), c.ntt(limit, 0);
            for (int i = 0; i < limit; i++) a[i] = 1ll * a[i] * c[i] % P;
            a.ntt(limit, 1);
            b.set(d - 1);
            for (int i = 0; i < d; i++) b[i] = 1ll * (a[i] + b[i]) * INV2 % P;
        }
        b.set(n - 1);
        return b;
    }
    Polynomial derivative() {
        Polynomial b(len - 1);
        for (int i = 1; i <= len; i++) b[i - 1] = 1ll * a[i] * i % P;
        return b;
    }
    Polynomial integral() {
        Polynomial b(len + 1);
        for (int i = 0; i <= len; i++) b[i + 1] = 1ll * a[i] * qpow(i + 1, P - 2) % P;
        return b;
    }
    Polynomial ln(int n) {
        return (derivative() * inv(n)).set(n - 1).integral().set(n - 1);
    }
    Polynomial exp(int n) {
        Polynomial b;
        b[0] = 1;
        for (int d = 1; d < (n << 1); d <<= 1) {
            Polynomial a = *this - b.ln(d) + 1;
            a.set(d - 1);
            int limit = d << 1;
            calcRev(limit);
            a.ntt(limit, 0), b.ntt(limit, 0);
            for (int i = 0; i < limit; i++) {
                b[i] = 1ll * a[i] * b[i] % P;
            }
            b.ntt(limit, 1);
            b.set(d - 1);
        }
        b.set(n - 1);
        return b;
    }
    Polynomial pow(int b, int n) {
        return (ln(n) * b).exp(n);
    }
    pair<Polynomial, Polynomial> operator/(Polynomial b) {
        int n = len, m = b.len;
        Polynomial ra = *this, rb = b, rc, d;
        ra.reverse(), rb.reverse();
        rc = ra * rb.inv(n - m + 1);
        rc.set(n - m);
        rc.reverse();
        d = *this - b * rc;
        d.set(m - 1);
        return make_pair(rc, d);
    }
};