我们将 $A,B$ 的每一位看成两个多项式的系数，那么对 $A,B$ 进行多项式乘法就可以得到乘积对应的多项式的各次项系数。最后进位即可求得答案。

下面探讨两个多项式进行乘法的问题。

列竖式

回顾我们小学二年级是如何做多项式乘法的？一般的做法应该是枚举两个多项式的每一项相乘，并将结果相加。这样一来，如果两个多项式长度分别为 $n,m$，一次乘法需要做 $\mathcal{O}(nm)$ 次运算。当 $n,m$ 特别大的时候，这是无法接受的。

那有什么更高效的做法吗？有的兄弟，有的！

前置 1

复数的三角表示

对于任意一个复数 $z$，我们可以将其表示成 $z=r(\cos\theta+i\sin \theta)$ 的形式，其中 $r=|z|$，$\theta$ 为 $z$ 在的复平面上对应的向量与实轴的夹角。

基于这种表示方法，我们可以发现对于两个复数 $z_1=r_1(\cos\theta_1+i\sin\theta_1),z_2=r_2(\cos\theta_2+i\sin\theta_2)$，有：

$$\begin{aligned}z_1z_2&=r_1r_2[\cos\theta_1\cos\theta_2-\sin\theta_1\sin\theta_2+i(\sin\theta_1\cos\theta_2+\cos\theta_1\sin\theta_2)]\\&=r_1r_2[\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2)]\end{aligned} $$

我们将其概括为：模长相乘、辐角相加。

可以这样实现一个复数类：

struct FS {
    db x, y; FS(db a = 0, db b = 0) { x = a; y = b; }
    FS operator+(const FS &o) const { return FS(x + o.x, y + o.y); }
    FS operator-(const FS &o) const { return FS(x - o.x, y - o.y); }
    FS operator*(const FS &o) const {
        return FS(x * o.x - y * o.y, x * o.y + y * o.x);
    }
};

单位根

方程 $x^n=1(n\in \mathbb{N+})$ 的全体复数解被称为 $n$ 次单位根。由代数基本定理可知，这样的解只有 $n$ 个。根据模长相乘、辐角相加，可以发现满足 $z=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}(k\in[0,n)\cap \mathbb{Z})$ 的复数 $z$ 均是 $n$ 次单位根且辐角互不相同。那么我们就这样找到了 $n$ 个本质不同的单位根。

注意到这些单位根将单位圆 $n$ 等分。我们通常记 $\omega_n$ 表示第一个单位根即 $\omega_n=\cos\dfrac{2\pi}{n}+i\sin\dfrac{2\pi}{n}$。则 $n$ 个单位根可以记为 $\omega_n^0,\dots,\omega_n^{n-1}$，因为 $\omega_{n}^k=\left(\cos\dfrac{2\pi}{n}+i\sin\dfrac{2\pi}{n}\right)^k=\cos\dfrac{2k\pi}{n}+i\sin\dfrac{2k\pi}{n}$。

单位根有如下性质：

• $\omega_n^n=1$。由定义可知。
• $\omega_n^k=\omega_{2n}^{2k}$。代入 $\omega_{n}^k$ 和 $\omega_{2n}^{2k}$ 的表达式后可以发现三角函数中的角度是相同的。
• $\omega_{2n}^{k+n}=-\omega_{2n}^k$。代入两者表达式后可以发现三角函数中的角度相差 $\pi$，由诱导公式可知实部和虚部都分别取相反数。

多项式的点值表示

$n$ 次多项式的一般形式为 $A(x)=\sum\limits_{i=0}^na_ix^i$。而如果我们取不同的 $n+1$ 个实数 $x_1,\dots,x_{n+1}$，且知道将她们代入多项式后多项式的值 $y_1,\dots,y_{n+1}$，那么可以唯一确定一个不超过 $n$ 次的多项式。此处人为规定保证存在这样的多项式。此时称 $(x_i,y_i)$ 为点值，称 $\{(x_1,y_1),\dots,(x_{n+1},y_{n+1})\}$ 为多项式的一个点值表示。一个多项式可以有多种点值表示，但是一个点值表示对应的多项式是唯一的。

考虑反证。假设存在两个多项式 $A,B$ 都满足这样的点值表示，那么进行移项可以得到一个不超过 $n$ 次的多项式 $C$。而每个 $x_i$ 都是这个多项式的零点，考虑展开其因式分解形式后一定会出现不低于 $n+1$ 次项，矛盾。

因此，$n+1$ 个点值可以确定至多一个不超过 $n$ 次的多项式。

快速傅里叶变换（FFT）

分治法

正片开始。考虑两个最高次数分别为 $n,m$ 的多项式 $A,B$。如果用系数表示，则难以合并积的系数。但是用点值表示时容易合并积的点值。因此考虑求出积的一个点值表示。

我们将 $A,B$ 两个多项式的长度补成 $2^k$，其中 $2^k$ 是最小的不少于 $n+m+1$ 的 $2$ 的幂。接下来令 $N=2^k$。考虑求出 $N$ 个点值，问题是选取哪 $N$ 个点值呢？

答案是 $N$ 次单位根。考虑 $A(x)=\sum\limits_{i=0}^{N-1}a_ix^i,F(x)=\sum\limits_{i=0}^{\tfrac{N}{2}-1}a_{2i}x^i,G(x)=\sum\limits_{i=0}^{\frac{N}{2}-1}a_{2i+1}x^i$，则 $A(x)=F(x^2)+xG(x^2)$。就是把奇数次项和偶数次项系数分别提取出来构成新的多项式。那么对于一个 $N$ 次单位根 $\omega_N^k\left(k\in\left[0,\dfrac{N}{2}\right)\cap \mathbb{Z}\right)$，运用单位根的性质，有：

• $A\left(\omega_{N}^k\right)=F\left(\omega_N^{2k}\right)+\omega_{N}^kG\left(\omega_N^{2k}\right)=F\left(\omega_{\frac{N}{2}}^k\right)+\omega_{N}^{k}G\left(\omega_{\frac{N}{2}}^k\right)$。
• $A\left(\omega_N^{k+\frac{N}{2}}\right)=F\left(\omega_N^{2k+N}\right)+\omega_N^{k+\frac{N}{2}}G\left(\omega_N^{2k+N}\right)=F\left(\omega_{\frac{N}{2}}^k\right)-\omega_{N}^{k}G\left(\omega_{\frac{N}{2}}^k\right)$。

那么我们只需要求出 $F,G$ 两个多项式在 $\dfrac{N}{2}$ 次单位根的点值表示即可。注意到 $F,G$ 的长度均为 $\dfrac{N}{2}$，所以这是相同两个子问题，直接递归计算即可。

每层递归子问题规模除以 $2$，所以 $T(n)=2T\left(\dfrac{n}{2}\right)+\mathcal{O}(n)$。根据主定理，时间复杂度为 $\mathcal{O}\left(N\log N\right)$。

void FFT(FS *a, int n) {
    if (n == 1) return; int m = n >> 1; FS f[m], g[m];
    for (int i = 0; i < n; i += 2) f[i >> 1] = a[i], g[i >> 1] = a[i | 1];
    FFT(f, m); FFT(g, m); FS w(1, 0), wn(cos(pi / m), sin(pi / m));
    for (int i = 0; i < m; ++i, w = w * wn)
        a[i] = f[i] + w * g[i], a[i + m] = f[i] - w * g[i];
}

代码中是传入多项式的系数然后将点值直接覆写在对应的系数上。

倍增法

引入

分治法的每层递归都要额外开销 $\mathcal{O}(N)$ 的空间。空间复杂度也为 $\mathcal{O}(N\log N)$。有没有什么优化的办法呢？

有的兄弟，有的！

位逆序置换

考虑用类似线段树的结构去刻画这个分治树的形态。令 $N=2^k$，定义序列 $p_j(i)$：

• $j=k$ 时，$p_j(i)=i$。
• $j\in [0,k)\cap\mathbb{Z}$ 时，记 $i=x2^{j+1}+y,y\in[0,2^{j+1})\cap\mathbb{Z}$。
  • $y\in\left[0,2^j\right)\cap\mathbb{Z}$ 时，$p_j(i)=p_{j+1}\left(x2^{j+1}+2y\right)$，即 $p_{j+1}$ 在 $\left[x2^{j+1},(x+1)2^{j+1}\right)$ 中的第 $y$ 个偶数位置上的数（0-indexed）。
  • $y\in \left[2^{j},2^{j+1}\right)$ 时，$p_j(i)=p_{j+1}\left(x2^{j+1}+2(y-2^j)+1\right)$，即 $p_{j+1}$ 在 $\left[x2^{j+1},(x+1)2^{j+1}\right)$ 中的第 $y-2^j$ 个奇数位置上的数（0-indexed）。

容易发现 $p_j$ 都是 $0\sim N - 1$ 的排列。于是记录其逆排列为 $q_j$。

记原来多项式各次项系数为 $f_i$，以及多项式 $F_{i,l}(x)=\sum\limits_{j=0}^{2^i-1}f_{p_{i}(l+j)}x^j$，此处 $l=k2^i$。

那么：

$$F_{i-1,l}(x)=\sum\limits_{j=0}^{2^{i-1}-1}f_{p_{i-1}(l+j)}x^j $$$$F_{i-1,l+2^{i-1}}=\sum\limits_{j=0}^{2^{i-1}-1}f_{p_{i-1}(l+2^{i-1}+j)}x^j $$

根据上述构造方式，$p_{i-1}(l+j)=p_i(l+2j),p_{i-1}(l+2^{i-1}+j)=p_i(l+2j+1)$。相当于把 $F_{i,l}(x)$ 的奇偶次项系数分别提取，所以有：

$$F_{i,l}(x)=F_{i-1,l}\left(x^2\right)+xF_{i-1,l+2^{i-1}}\left(x^2\right) $$

根据 FFT 中的推导，$F_{i,l}(x)$ 在 $2^i$ 个单位根处的点值可以由 $F_{i-1,l}(x)$ 和 $F_{i-1,l+2^{i-1}}$ 在 $2^{i-1}$ 个单位根处的点值计算得出。在分治法中，我们递归两个多项式求解子问题再推出原问题的答案。于是在这里我们考虑由子问题答案直接推出原问题答案。

通过观察，我们发现对于 $q_0$ 满足：把 $i$ 和 $q_0(i)$ 都表示成 $k$ 位二进制数 $A=\sum\limits_{t=0}^{k-1}a_t2^t,B=\sum\limits_{t=0}^{k-1}b_t2^t$，有 $a_t=b_{k-1-t}$。

证明：

考虑证明一个更强的命题：

对于任意 $j\in [0,k]$，记 $q_j(i)$ 的 $k$ 二进制表示为 $A_j(i)=\sum\limits_{t=0}^{k-1}a_j(t)2^t$，则有：

• $\forall \,t\in[0,k-j)\cap\mathbb{Z},b_t=a_j(k-1-t)$。
• $\forall\, t\in[k-j,k)\cap\mathbb{Z},b_t=a_j(t+j-k)$。

原命题即该命题在 $j=0$ 时的情况。

考虑使用数学归纳法。

• 当 $j=k$ 时，显然成立。

• 假设当 $j=u$ 时成立，则 $j=u-1$ 时：

  记 $q_u(i)=X2^u+Y,Y\in \left[0,2^u\right)\cap\mathbb{Z}$，以及 $q_u(i)\bmod 2=Y\bmod 2 = v,v\in\{0,1\}$。

  根据构造方式，我们发现 $p_{u-1}\left(X2^u+v2^{u-1}+\dfrac{Y-v}{2}\right)=p_u(q_u(i))=i$，因此 $q_{u-1}(i)=X2^u+v2^{u-1}+\dfrac{Y-v}{2}$。

  熟悉位运算的话可能已经看出来了，就是把原来 $2,\dots,2^{u-1}$ 这些二进制位右移一位，再把原来最低位上的数放到第 $2^{u-1}$ 这个位置上。但是为了不失严谨性，考虑继续设 $X=\sum\limits_{t=0}^{k-u}x_t2^t,Y=\sum\limits_{t=0}^{u-1}y_t2^t$，则 $v=y_0$。

  那么有：

  $$q_u(i)=\sum\limits_{t=0}^{u-1}y_t2^t+\sum\limits_{t=u}^kx_{t-u}2^t $$

  和：

  $$\begin{aligned}q_{u-1}(i)&=\sum\limits_{t=1}^{u-1}y_t2^{t-1}+y_02^{u-1}+\sum\limits_{t=u}^kx_{t-u}2^t\\&=\sum\limits_{t=0}^{u-2}y_{t+1}2^t+y_02^{u-1}+\sum\limits_{t=u}^kx_{t-u}2^t\end{aligned} $$

  那么记 $q_u(i)$ 和 $q_{u-1}(i)$ 的二进制表示分别为 $\sum\limits_{t=0}^{k-1}\alpha_t2^t$ 和 $\sum\limits_{t=0}^{k-1}\beta_t2^t$，与上面两式对应得：

  $$\alpha_t=\begin{cases}y_t,t\in[0,u)\cap\mathbb{Z}\\x_{t-u},t\in [u,k)\cap\mathbb{Z}\end{cases} $$$$\beta_t=\begin{cases}y_{t+1},t\in[0,u-1)\cap\mathbb{Z}\\y_0,t=u-1\\x_{t-u},t\in [u,k)\cap\mathbb{Z}\end{cases} $$

  那么：

  • 当 $t\in[0,k-u+1)\cap\mathbb{Z}$ 时：
    • 当 $t\in[0,k-u)\cap\mathbb{Z}$ 时，$k-1-t\in [u,k)\cap\mathbb{Z}$。注意到此时 $\alpha_{k-1-t}=\beta_{k-1-t}$，且由归纳假设知 $\alpha_{k-1-t}=b_t$，因此 $\beta_{k-1-t}=b_t$。
    • 当 $t=k-u$ 时，$k-1-t=u-1$。此时 $\beta_{k-1-t}=y_0=\alpha_0$。由归纳假设知 $\alpha_0=b_{k-u}=b_t$，所以 $\beta_{k-1-t}=b_t$。
  • 当 $t\in[k-u+1,k)\cap\mathbb{Z}$ 时，$t+u-1-k\in [0,u-1)\cap\mathbb{Z}$。可知此时 $\beta_{t+u-1-k}=\alpha_{t+u-k}$。由归纳假设可知 $\alpha_{t+u-k}=b_t$，所以 $\beta_{t+u-1-k}=b_t$。

  和归纳假设的形式一模一样。故证毕。

考虑求出每个位置 $i$ 的逆序 $R_i$。$\mathcal{O}(n\log n)$ 是容易的。不过可以线性解决。

考虑 $i=\sum\limits_{j=0}^{k-1}a_j2^j$，那么 $\left\lfloor\dfrac{i}{2}\right\rfloor=0\times 2^{k-1}+\sum\limits_{j=0}^{k-2}a_{j+1}2^j$。

此时 $R_i=\sum\limits_{j=0}^{k-1}a_{k-1-j}2^j,R_{\lfloor\frac{i}{2}\rfloor}=\sum\limits_{j=1}^{k-1}a_{k-j}2^j+0$。注意到：

$$R_i=a_{0}2^{k-1}+\sum\limits_{j=0}^{k-2}a_{k-1-j}2^j=a_02^{k-1}+\sum\limits_{j=1}^{k-1}a_{k-j}2^{j-1}=a_02^{k-1}+\dfrac{R_{\lfloor\frac{i}{2}\rfloor}}{2} $$

那么可以直接 $\mathcal{O}(n)$ 递推出所有 $R_i$。

蝶形运算优化

通过位逆序置换，我们已经得到所有 $F_{0,i}(x)$ 的答案了。而我们的目标是求出 $F_{k,0}(x)$ 的答案。

考虑将点值直接覆写在 $f$ 数组上。更具体地说，进行 $m$ 轮覆写，第 $i$ 轮覆写中，对于 $F_{i,l}(x)$，我们将 $f_{l+t}(t\in [0,2^i)\cap\mathbb{Z})$ 覆写为 $F_{i,l}\left(\omega_{2^i}^t\right)$。

类似数学归纳法。位逆序置换相当于覆写 $m=0$ 时的答案。假设我们已经覆写好了第 $i-1$ 轮的答案。那么在第 $i$ 轮时，对于 $t\in\left[0,2^{i-1}\right)\cap\mathbb{Z}$，根据 FFT 的结论：

• $F_{i,l}\left(\omega_{2^i}^t\right)=F_{i-1,l}\left(\omega_{2^{i-1}}^t\right)+\omega_{2^i}^tF_{i-1,l+2^{i-1}}\left(\omega_{2^{i-1}}^t\right)$。
• $F_{i,l}\left(\omega_{2^i}^{t+2^{i-1}}\right)=F_{i-1,l}\left(\omega_{2^{i-1}}^t\right)-\omega_{2^i}^tF_{i-1,l+2^{i-1}}\left(\omega_{2^{i-1}}^t\right)$。

根据覆写的结果，$F_{i,l}\left(\omega_{2^i}^t\right)=f_{l+t}+\omega_{2^i}^tf_{l+2^{i-1}+t},F_{i,l}\left(\omega_{2^i}^{t+2^{i-1}}\right)=f_{l+t}-\omega_{2^i}^tf_{l+2^{i-1}+t}$。

那么令 $u=f_{l+t},v=\omega_{2^i}^tf_{l+2^{i-1}+t}$。实现的时候枚举每个长度为 $2^i$ 的段，再从小往大枚举 $t$ 进行覆写，那么每个 $t$ 时只有比 $t$ 小的位置被覆写，其她位置仍然是上一轮的覆写结果。那么记录下来的 $u,v$ 就是我们在这一轮覆写中需要用到的值。令 $f_{l+t}\leftarrow u+v,f_{l+2^{i-1}+t}\leftarrow u-v$ 即可完成第 $i$ 轮的覆写操作。

结合以上两个优化，可以得到非递归版 FFT 的代码：

void init(int n) {
	p[0] = 0;
	for (int i = 1; i < n; ++i) {
		p[i] = p[i >> 1] >> 1; if (i & 1) p[i] |= n >> 1;
	}
}
void FFT(FS *a, int n) {
    for (int i = 0; i < n; ++i) if (i < p[i]) swap(a[i], a[p[i]]);
	for (int i = 2; i <= n; i <<= 1) {
		FS g(cos(2 * pi / i), sin(2 * pi / i));
		for (int j = 0; j < n; j += i) {
			FS w(1, 0), u, v;
			for (int k = j; k < j + (i >> 1); ++k)
				u = a[k], v = w * a[k + (i >> 1)],
				a[k] = u + v, a[k + (i >> 1)] = u - v, w = w * g;
		}
	}
}

时间复杂度不变，空间复杂度优化为线性，常数也变小。

快速傅里叶逆变换

我们成功地得到了多项式的点值表示，那么如何将这个点值表示还原成系数表示呢？

记 $A\left(\omega_N^k\right)=y_k$。我们考虑构造一个新的多项式 $B(x)=\sum\limits_{i=0}^{N-1}y_ix^i$。

考察这个多项式在 $\omega_N^{-k}$ 处的点值：

$$\begin{aligned}B\left(\omega_N^{-k}\right)&=\sum\limits_{i=0}^{N-1}\omega_N^{-ki}\sum\limits_{j=0}^{N-1}a_j\omega_N^{ij}\\&=\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1}a_j\omega_N^{i(j-k)}\\&=\sum\limits_{j=0}^{N-1}\sum\limits_{i=0}^{N-1}a_j\left(\omega_N^{j-k}\right)^i\\&=\sum\limits_{j=0}^{N-1}a_j\sum\limits_{i=0}^{N-1}\left(\omega_N^{j-k}\right)^i\end{aligned} $$

记 $C(k)=\sum\limits_{i=0}^{N-1}\left(\omega_N^k\right)^i$，其中 $k\in[0,N)\cap\mathbb{Z}$。

• 当 $k=0$ 时，$C(k)=1$。
• 当 $k\ne 0$ 时，运用等比数列求和公式可知：$$C(k)=\dfrac{1-\omega_N^N}{1-\omega_N}=0$$

代回原式，$B\left(\omega_N^{-k}\right)=\sum\limits_{j=0}^{N-1}a_jC(j-k)=Na_k$。因为只有在 $j=k$ 时 $C(j-k)$ 才不为 $0$。

也就是说 $Na_k$ 是 $B(x)$ 在 $\omega_N^{-k}$ 处的值。注意到 $\omega_N^k=\omega_{N}^{N-k}$，那么单位根的指数就全部在 $[0,N)$ 之间了。我们只需要对 $B(x)$ 求 $N$ 次单位根的点值即可得到 $A(x)$ 的各次项系数的 $N$ 倍，那么就能求出 $A(x)$ 了。

【模板】多项式乘法（FFT）

运用快速傅里叶（逆）变换，我们可以 $\mathcal{O}(N\log N)$ 时间求出两个多项式 $A(x),B(x)$ 的乘积。流程如下：

• 将 $A(x),B(x)$ 的长度补成 $2$ 的幂次，且不小于乘积的长度。
• 对 $A(x),B(x)$ 求 $N$ 次单位根处的值。然后把对应点值乘起来就能得到 $A(x)B(x)$ 在 $N$ 次单位根处的值。
• 再以这些点值构造一个新的多项式，求出在 $N$ 次单位根处的值，即可还原出 $A(x)B(x)$ 的系数的 $N$ 倍。再将这些值除以 $N$ 即可求出答案。

由于实数运算会有精度损失，因此输出的时候要加上一个偏移量再取整。

递归版代码（2.31s 141.02 MB）：

#include <bits/stdc++.h>
#define db double
using namespace std; const db pi = acos(-1); const int N = 2.5e6;
int n, m, L;
struct FS {
    db x, y; FS(db a = 0, db b = 0) { x = a; y = b; }
    FS operator+(const FS &o) const { return FS(x + o.x, y + o.y); }
    FS operator-(const FS &o) const { return FS(x - o.x, y - o.y); }
    FS operator*(const FS &o) const {
        return FS(x * o.x - y * o.y, x * o.y + y * o.x);
    }
} a[N], b[N];
void FFT(FS *a, int n) {
    if (n == 1) return; int m = n >> 1; FS f[m], g[m];
    for (int i = 0; i < n; i += 2) f[i >> 1] = a[i], g[i >> 1] = a[i | 1];
    FFT(f, m); FFT(g, m); FS w(1, 0), wn(cos(pi / m), sin(pi / m));
    for (int i = 0; i < m; ++i, w = w * wn)
        a[i] = f[i] + w * g[i], a[i + m] = f[i] - w * g[i];
}
signed main() {
    cin.tie(0), cout.tie(0), ios::sync_with_stdio(0); cin >> n >> m;
    for (int i = 0; i <= n; ++i) cin >> a[i].x;
    for (int i = 0; i <= m; ++i) cin >> b[i].x;
    for (L = 1; L <= n + m; L <<= 1); FFT(a, L); FFT(b, L);
    for (int i = 0; i < L; ++i) a[i] = a[i] * b[i]; FFT(a, L);
    for (int i = 0; i <= n + m; ++i)
        cout << int((i ? a[L - i].x : a[i].x) / L + 0.5) << ' ';
    return 0;
}

非递归版代码（1.82s 84.91MB）：

#include <bits/stdc++.h>
#define db double
using namespace std; const db pi = acos(-1); const int N = 2.5e6;
int n, m, L, p[N];
struct FS {
    db x, y; FS(db a = 0, db b = 0) { x = a; y = b; }
    FS operator+(const FS &o) const { return FS(x + o.x, y + o.y); }
    FS operator-(const FS &o) const { return FS(x - o.x, y - o.y); }
    FS operator*(const FS &o) const {
        return FS(x * o.x - y * o.y, x * o.y + y * o.x);
    }
} a[N], b[N];
void init(int n) {
	p[0] = 0;
	for (int i = 1; i < n; ++i) {
		p[i] = p[i >> 1] >> 1; if (i & 1) p[i] |= n >> 1;
	}
}
void FFT(FS *a, int n) {
    for (int i = 0; i < n; ++i) if (i < p[i]) swap(a[i], a[p[i]]);
	for (int i = 2; i <= n; i <<= 1) {
		FS g(cos(2 * pi / i), sin(2 * pi / i));
		for (int j = 0; j < n; j += i) {
			FS w(1, 0), u, v;
			for (int k = j; k < j + (i >> 1); ++k)
				u = a[k], v = w * a[k + (i >> 1)],
				a[k] = u + v, a[k + (i >> 1)] = u - v, w = w * g;
		}
	}
}
signed main() {
    cin.tie(0), cout.tie(0), ios::sync_with_stdio(0); cin >> n >> m;
    for (int i = 0; i <= n; ++i) cin >> a[i].x;
    for (int i = 0; i <= m; ++i) cin >> b[i].x;
    for (L = 1; L <= n + m; L <<= 1); init(L); FFT(a, L); FFT(b, L);
    for (int i = 0; i < L; ++i) a[i] = a[i] * b[i]; FFT(a, L);
    for (int i = 0; i <= n + m; ++i)
        cout << int((i ? a[L - i].x : a[i].x) / L + 0.5) << ' ';
    return 0;
}

卷积

FFT 究竟干了什么

对于两个多项式 $A(x)=\sum\limits_{i=0}^na_ix^i,B(x)=\sum\limits_{i=0}^mb_ix^i$，我们求出了第 $A(x)B(x)$ 的各项系数。不妨认为 $n>m$，并将 $A(x),B(x)$ 用 $0$ 补至 $2n$ 次项系数。

考察 $A(x)B(x)$ 的第 $i$ 项系数究竟是什么：

$$\begin{aligned}A(x)B(x)&=\sum\limits_{i=0}^na_ix^i\sum\limits_{j=0}^nb_jx^j\\&=\sum\limits_{i=0}^n\sum\limits_{j=0}^na_ib_jx^{i+j}\\&=\sum\limits_{i=0}^n\sum\limits_{j=0}^n\sum\limits_{k=0}^{2n}[i+j=k]a_ib_jx^k\\&=\sum\limits_{k=0}^{2n}\sum\limits_{i=0}^n\sum\limits_{j=0}^n[i+j=k]a_ib_jx^k\\&=\sum\limits_{k=0}^{2n}\sum\limits_{i=0}^ka_ib_{k-i}x^k\\&=\sum\limits_{k=0}^{2n}\left(\sum\limits_{i=0}^ka_ib_{k-i}\right)x^k\end{aligned} $$

牢记这个卷积形式：

$$\color{red}{\sum\limits_{i=0}^na_ib_{n-i}}$$

特征是下标和为定值。这种类型的式子可以通过构造多项式乘法进行计算。下标差为定值时，可以通过构造一个数组的对称数组，从而转化为下标和为定值的情况。

习题 P3338 [ZJOI2014] 力

给出 $[q_1,\dots,q_n]$，对于每一个 $j$，求：

$$E_j=\sum\limits_{i=1}^{j-1}\dfrac{q_j}{(i-j)^2}-\sum\limits_{i=j+1}^n\dfrac{q_j}{(i-j)^2} $$

$n\le 10^5$。精度误差不超过 $10^{-2}$。

先推式子，令 $k=|i-j|$：

$$\begin{aligned}E_j&=\sum\limits_{k=1}^{j-1}\dfrac{q_{j-k}}{k^2}-\sum\limits_{k=1}^{n-j}\dfrac{q_{j+k}}{k^2}\end{aligned} $$

令 $g_i=\dfrac{1}{i^2}$：

$$\begin{aligned}E_j&=\sum\limits_{k=1}^{j-1}{q_{j-k}}g_k-\sum\limits_{k=1}^{n-j}q_{j+k}g_k\end{aligned} $$

前半部分已经是卷积形式，可以看成是一次项到 $n$ 次项分别为 $q_i,g_i$ 的两个多项式卷积后的第 $j$ 项系数。

令 $f_i=q_{n-i+1}$：

$$\begin{aligned}\sum\limits_{k=1}^{n-j}q_{j+k}g_k=\sum\limits_{k=1}^{n-j}f_{n-j+1-k}g_k\end{aligned} $$

那么也被我们化成了 $f_i,g_i$ 为各次项系数的两个多项式卷积后的第 $n-j+1$ 次项系数。

FFT 计算即可。时空复杂度 $\mathcal{O}(n\log n)$。

#include <bits/stdc++.h>
using namespace std; typedef long double db;
const int N = 1 << 18; const db pi = acos(-1); int n, L = 1; db q[N];
struct FS {
	db x, y; FS(db a = 0, db b = 0) { x = a; y = b; }
	FS operator+(const FS &o) const { return FS(x + o.x, y + o.y); }
	FS operator-(const FS &o) const { return FS(x - o.x, y - o.y); }
	FS operator*(const FS &o) const {
		return FS(x * o.x - y * o.y, x * o.y + y * o.x);
	}
} a[N], b[N], c[N], d[N];
void FFT(FS *a, int n) {
	if (n == 1) return; int m = n >> 1; FS f[m], g[m];
	for (int i = 0; i < m; ++i) f[i] = a[i << 1], g[i] = a[i << 1 | 1];
	FFT(f, m); FFT(g, m); FS w(1, 0), k(cos(pi / m), sin(pi / m));
	for (int i = 0; i < m; ++i, w = w * k)
		a[i] = f[i] + w * g[i], a[i + m] = f[i] - w * g[i];
}
void conv(FS *a, FS *b, FS *c, int n) {
	FFT(a, n); FFT(b, n); for (int i = 0; i < n; ++i) a[i] = a[i] * b[i];
	FFT(a, n); c[0].x = a[0].x / n;
	for (int i = 1; i < n; ++i) c[i].x = a[n - i].x / n;
}
int main() {
	cin.tie(0); cout.tie(0); ios::sync_with_stdio(0); cin >> n;
	for (int i = 1; i <= n; ++i)
		cin >> q[i], a[i] = FS(q[i], 0), b[i] = FS(db(1.0 / i / i), 0);
	while (L <= (n << 1)) L <<= 1; conv(a, b, c, L); a[0] = b[0] = FS();
	for (int i = 1; i <= n; ++i)
		a[i] = FS(q[n - i + 1], 0), b[i] = FS(db(1.0 / i / i), 0);
	for (int i = n + 1; i < L; ++i) a[i] = b[i] = FS(); conv(a, b, d, L);
	for (int i = 1; i <= n; ++i)
		cout << fixed << setprecision(3) << c[i].x - d[n - i + 1].x << '\n';
	return 0;
}

前置 2

阶

定义

由欧拉定理可知，若 $\gcd(a,m)=1$，则 $a^{\varphi(m)}\bmod m=1$。

因此满足 $a^n\bmod m=1$ 的最小正整数 $n$ 存在，这个 $n$ 称作 $a$ 模 $m$ 的阶，记为 $\delta_m(a)$。

性质 1

$\forall\, i,j\in[1,\delta_m(a)]\cap\mathbb{Z},a^i\bmod m\ne a^j\bmod m$。

证明：

不妨认为 $i<j$。考虑反证，若 $\exist\,i,j\in [1,\delta_m(a)]\cap\mathbb{Z},a^i\bmod m= a^j\bmod m$，那么 $a^i(a^{j-i}-1)\bmod m =0$，即 $a^i(a^{j-i}-1)$ 是 $m$ 的倍数。

由裴蜀定理得 $\gcd(a^i,m)=1$，再次由其得 $a^{j-i}-1$ 是 $m$ 的倍数，故 $a^{j-i}\bmod m=1$，那么存在一个更小的正整数 $k=j-i$ 使得 $a^k\bmod m=1$，这与阶的定义矛盾。

故原命题得证。

性质 2

若 $a^n\bmod m =1$，则 $n\bmod \delta_m(a)=0$。

证明：

设 $n=k\delta_m(a)+r,r\in[0,\delta_m(a))\cap \mathbb{Z}$。

考虑反证，若 $r>0$，则：

$$a^r\bmod m=a^r\left(a^{\delta_m(a)}\right)^k\bmod m=a^n\bmod m = 1 $$

所以存在更小的正整数 $r$ 使得 $a^r\bmod m=1$，这与阶的定义矛盾。

故原命题得证。

原根

定义

若 $\gcd(g,m)=1$ 且 $\delta_m(g)=\varphi(m)$，则称 $g$ 是模 $m$ 的原根。

判定定理

设 $m\ge 3,\gcd(g,m)=1$，则 $g$ 是模 $m$ 的原根的充要条件是对于 $\varphi(m)$ 的每个质因数 $p$，都有 $g^{\frac{\varphi(m)}{p}}\bmod m\ne 1$。

证明：

必要性显然，下证充分性。

考虑反证法，假设这样的 $g$ 不是模 $m$ 的原根。

由于 $g^{\varphi(m)}\bmod m=1$，由性质 2 得 $\varphi(m)\bmod \delta_m(g)=0$。且 $g$ 不是模 $m$ 的原根，因此 $\delta_m(g)<\varphi(m)$。

记 $k=\dfrac{\varphi(m)}{\delta_m(g)}$，可知 $k\ge 2$，其存在质因数 $p$。记 $k=pq$，则 $\varphi(m)=pq\delta_m(g)$。此时有：

$$g^{\frac{\varphi(m)}{p}}\bmod m =\left(g^{\delta_m(g)}\right)^q\bmod m=1 $$

故存在 $\varphi(m)$ 的质因数 $p$ 使得 $g^{\frac{\varphi(m)}{p}}\bmod m = 1$，与条件矛盾。

故原命题得证。

快速数论变换（NTT）

用原根代替单位根

在 FFT 中，我们需要用到 double 类型，会导致计算结果的精确度下降。那么有没有解决方法呢？

有的兄弟，有的！

考虑一个质数 $p=qN+1$，其中 $N=2^m$。那么模 $p$ 的原根 $g$ 满足 $g^{qN}\bmod p=1$，且由性质 1 可知 $g^i(i\in[0,qN]\cap\mathbb{Z})$ 在模 $m$ 意义下两两不同。

对于 $n=2^m$，定义 $n$ 阶模 $p$ 单位根 $w_n$（名字和记号是我瞎起的，为了与单位根 $\omega_n$ 区分）为 $g^{\frac{qN}{n}}$。$w_n$ 有如下性质：

• $w_n^n\bmod p = 1$。由定义可知。
• $w_n^k=w_{2n}^{2k}$。因为 $w_{2n}^{2k}=g^{\frac{2kqN}{2n}}=g^{\frac{kqN}{n}}=w_n^k$。
• $w_{2n}^{n}\bmod p= p-1$。因为 $\left(w_{2n}^n\right)^2\bmod p=w_{2n}^{2n}\bmod p=1$，所以 $(w_{2n}^n+1)(w_{2n}^n-1)\bmod p = 0$。根据性质 1 可知 $w_{2n}^n\bmod p\ne 1$，所以 $(w_{2n}^n-1)\bmod p\ne 0$，因此 $\gcd(w_{2n}^n-1,p)=1$，由裴蜀定理得 $(w_{2n}^n+1)\bmod p = 0$，即 $w_{2n}^n\bmod p = p-1$。
• $(w_{2n}^k+w_{2n}^{k+n})\bmod p = 0$。容易由上一条推出。

我们发现 $w_n$ 与 $\omega_n$ 有相似的性质！

那么考虑用 $w_n$ 替代 $\omega_n$ 进行模意义下的计算。此处推导只需要把 FFT 中的推导的等号全部改成模 $p$ 意义下相等即可。不再赘述。

值得注意的是，此时求出来的点值以及系数都是模 $p$ 意义下的。在数据范围较小（如模板题）时，求出来的系数就是真实值。

常见的质数和原根有：

• $p=167772161=5\times 2^{25}+1,g=3$。
• $p=469762049=7\times 2^{26}+1,g=3$。
• $p=754974721=45\times 2^{24}+1,g=11$。
• $\color{red} p=998244353=119\times 2^{23}+1,g=3$。熟悉吗？
• $p=1004535809=479\times 2^{21}+1,g=3$。

时间复杂度 $\mathcal{O}(n\log n)$，空间复杂度 $\mathcal{O}(n)$。1.38s 24.47MB。

#include <bits/stdc++.h>
template<class T> void read(T &x) {
    x = 0; T f = 1; char c = getchar();
    for (; c < 48 || c > 57; c = getchar()) if (c == 45) f = -1;
    for (; c > 47 && c < 58; c = getchar()) x = x * 10 + c - 48; x *= f;
}
template<class T> void write(T x) {
    if (x > 9) write(x / 10); putchar(x % 10 + 48);
}
template<class T> void print(T x, char ed = '\n') {
    if (x < 0) putchar('-'), x = -x; write(x), putchar(ed);
}
using namespace std; const int N = 1 << 21, M = 998244353, G = 3;
int n, m, a[N], b[N], L, p[N];
int qp(int x, int y) {
    int r = 1;
    for (; y; y >>= 1, x = 1ll * x * x % M) if (y & 1) r = 1ll * r * x % M;
    return r;
}
void init(int n) {
	p[0] = 0;
	for (int i = 1; i < n; ++i) {
		p[i] = p[i >> 1] >> 1; if (i & 1) p[i] |= n >> 1;
	}
}
void NTT(int *a, int n) {
	for (int i = 0; i < n; ++i) if (i < p[i]) swap(a[i], a[p[i]]);
	for (int i = 2, g; i <= n; i <<= 1) {
		g = qp(G, (M - 1) / i);
		for (int j = 0, w; j < n; j += i) {
			w = 1;
			for (int k = j, u, v; k < j + (i >> 1); ++k)
				u = a[k], v = 1ll * w * a[k + (i >> 1)] % M, a[k] = (u + v) % M,
				a[k + (i >> 1)] = (u - v + M) % M, w = 1ll * w * g % M;
		}
	}
}
signed main() {
    read(n); read(m); for (int i = 0; i <= n; ++i) read(a[i]);
    for (int i = 0; i <= m; ++i) read(b[i]);
    for (L = 1; L <= n + m; L <<= 1); init(L); NTT(a, L); NTT(b, L);
    for (int i = 0; i < L; ++i) a[i] = 1ll * a[i] * b[i] % M; NTT(a, L);
    for (int i = 0, j = qp(L, M - 2); i <= n + m; ++i)
        print(1ll * (i ? a[L - i] : a[0]) * j % M, ' ');
    return 0;
}

习题 P3723 [AH2017/HNOI2017] 礼物

• 给出 $[x_1,\dots,x_n],[y_1,\dots,y_n]$，对于任意整数 $c$ 和 $[y_1,\dots,y_n]$ 的循环同构 $[p_1,\dots,p_n]$，求下式最小值：

$$\sum\limits_{i=1}^n(x_i-p_i+c)^2$$

• $n\le 50000$，$x_i,y_i\le 100$。

将式子展开可以得到：

$$\sum\limits_{i=1}^n\left(x_i^2+p_i^2\right)+2\sum\limits_{i=1}^n(x_i-p_i)c+nc^2-2\sum\limits_{i=1}^nx_ip_i $$

注意到 $p$ 是 $y$ 的循环同构，因此一、二次方和相等。那么第一项是常数。第二、三项是一个关于 $c$ 的二次函数。此时变量 $c$ 和 $p$ 独立开来，分别求解。二次函数部分容易求出最小值。

然后我们希望求出 $\sum\limits_{i=1}^nx_ip_i$ 的最大值。考虑将 $y$ 倍长，那么如果 $p$ 是 $y$ 向左循环移位 $k(k\in[0,n))$ 位的结果，则 $p_i=y_{i+k}$。

现在要求 $\sum\limits_{i=1}^nx_iy_{i+k}$ 的最大值。考虑求出每个 $k$ 的答案。注意到若将 $k$ 看作定值，则此时下标差为定值 $k$，翻转 $y$ 序列得到 $[a_1,\dots,a_{2n}]$，则原式变成 $\sum\limits_{i=1}^{n}x_ia_{2n-k+1-i}$。

这样式子又变成了我们熟悉的卷积形式，NTT 计算即可。更具体地，令 $a_0=x_0=0$，则可以看成是 $F(t)=\sum\limits_{i=0}^nx_it^i,G(t)=\sum\limits_{i=0}^{2n}a_it^i$ 两个多项式卷积后的第 $2n-k+1$ 次项系数。其实这个系数和答案的求和范围有一些差距，但是由于 $x_i$ 只有在 $i\in [1,n]$ 时非零，因此两者相等。

由于数据范围不大，求出来模意义下的答案就是真实答案。

时间复杂度 $\mathcal{O}(n\log n)$，空间复杂度 $\mathcal{O}(n)$。

#include <bits/stdc++.h>
#define int long long
using namespace std; const int N = 1 << 19, M = 998244353, G = 3;
int n, m, a[N], b[N], c[N], p[N], L, ans, a1, sum, a0 = 0;
int qp(int x, int y) {
	int r = 1; for (; y; y >>= 1, x = x * x % M) if (y & 1) r = r * x % M; return r;
}
void NTT(int *a, int n, bool o = 0) {
	if (o)
		for (int i = 1; i < n; ++i) {
			p[i] = p[i >> 1] >> 1; if (i & 1) p[i] |= n >> 1;
		}
	for (int i = 0; i < n; ++i) if (i < p[i]) swap(a[i], a[p[i]]);
	for (int i = 2, g; i <= n; i <<= 1) {
		g = qp(G, (M - 1) / i);
		for (int j = 0, w; j < n; j += i) {
			w = 1;
			for (int k = j, u, v; k < j + (i >> 1); ++k, w = w * g % M)
				u = a[k], v = w * a[k + (i >> 1)] % M,
				a[k] = (u + v) % M, a[k + (i >> 1)] = (u - v + M) % M;
		}
	}
}
void conv(int *a, int n, int *b, int m, int *c, int &L) {
	int iv; for (L = 1; L <= n + m; L <<= 1); NTT(a, L, 1); NTT(b, L);
	for (int i = 0; i < L; ++i) a[i] = a[i] * b[i] % M; NTT(a, L);
	iv = qp(L, M - 2); c[0] = a[0] * iv % M;
	for (int i = 1; i < L; ++i) c[i] = a[L - i] * iv % M;
}
int Abs(int x) { return max(x, -x); }
int calc(int x) { return n * x * x + (a1 << 1) * x; }
signed main() {
	cin.tie(0); cout.tie(0); ios::sync_with_stdio(0); cin >> n >> m;
	for (int i = 1; i <= n; ++i)
		cin >> a[i], a[i + n] = a[i], sum += a[i] * a[i], a1 += a[i];
	for (int i = 1; i <= n; ++i) cin >> b[i], a1 -= b[i], sum += b[i] * b[i];
	reverse(a + 1, a + (n << 1 | 1)); conv(a, n << 1, b, n, c, L);
	for (int i = 0; i < n; ++i) a0 = max(a0, c[(n << 1) - i + 1]);
	ans = min({calc(Abs(a1) / n), calc(Abs(a1) / n + 1),
	           calc(-(Abs(a1) / n)), calc(-(Abs(a1) / n + 1))
	      });
	return cout << sum + ans - (a0 << 1), 0;
}