FFT 与 NTT 学习笔记

【概念】

点值：给定多项式 \(f(x)=a_0+a_1x+\cdots+a_{n-1}x^{n-1}\) 和 \(x_1\sim x_m\)，求 \(f(x_1),f(x_2),\dots,f(x_m)\)。（\(m=n\)）

求点值的算法一般是 \(O(n^2)\) 的，但对于特殊的多项式，可以做到更快。

---

插值：给定 \((x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\)，其中 \(x_0\sim x_{n-1}\) 互不相同。求一个不超过 \(n-1\) 次的多项式 \(f(x)\)，使得 \(f(x_i)=y_i\)。可以观察到 \(f(x)\) 是存在且唯一的。

求插值的算法一般是拉格朗日插值，令多项式 \(P_i(x)=\dfrac{\displaystyle\prod_{j=0,j\neq i}^{n-1}(x-x_j)}{\displaystyle\prod_{j=0,j\neq i}^{n-1}(x_i-x_j)}\)。

容易观察到 \(P_i(x_i)=1,P_i(x_j)=0\text{（j 不等于 i）}\)。

则 \(f(x)=\displaystyle\prod P_i(x_i)\cdot y_i\)。这证明了 \(f(x)\) 是存在的。

---

下证唯一性。若有两个多项式 \(f(x),g(x)\) 均满足要求。

考虑 \(r(x)=f(x)-g(x)\)，则 \(r(x)\) 在 \(x_0\sim x_{n-1}\) 处值为 \(0\)。

\(r(x)\) 是 \(n-1\) 次多项式，却有 \(n\) 个不相同的根，所以 \(r(x)=0\)，所以 \(f(x)=g(x)\)。

【单位根的性质】

\(n\) 次单位根是所有满足 \(z^n=1\) 的 \(z\) 们。

考虑复数，记 \(n\) 次单位根们为 \(\varepsilon_n\)。（即在复单位圆上从 \(0i+1\) 逆时针转 \(\dfrac{360°}{n}\) 得到的复数）

因为单位根的定义，所以 \(\varepsilon_n^{0\sim n-1}\) 都满足 \(z^n=1\)，所以他们刚好构成 \(z^n=1\) 的所有解。

还有一个性质：\(\varepsilon_{n}^{k}=e^{i\cdot \frac{2\pi}{n}\cdot k}=\cos \dfrac{2\pi k}{n}+i\cdot \sin \dfrac{2\pi k}{n}\)。

复数的折半定理：\(\varepsilon_{2n}^{2i}=\varepsilon_n^{i}\)。

【FFT：快速傅里叶变换】

考虑这样一个问题：给定 \(f(x),g(x)\)，求 \(h(x)=f(x)\cdot g(x)\)。

多项式的乘积，其实就是系数的卷积。所以 FFT 解决的是卷积的问题。

令 \(n=\deg f(x) + \deg g(x) + 1\)，显然 \(\deg h(x)<n\)。

如果直接循环做，是 \(O(n^2)\) 的；但是 FFT 能做到 \(O(n\log n)\)。

---

为了确定 \(h(x)\)，只需要 \(h\) 在 \(n\) 个点 \(x_0\sim x_{n-1}\) 上的值即可。

基本思路是：① 求出 \(f(x_0)\sim f(x_{n-1})\)。 ② 求出 \(g(x_0)\sim g(x_{n-1})\)。 ③ \(h(x_i)=f(x_i)\cdot g(x_i)\) ④ 插值还原 \(h\)。

拉格朗日插值也是 \(O(n^2)\) 的，所以我们要用更快的插值方式，后面会说到。

这里我们取 \(x_i=\varepsilon_{n}^{i}\)。

（当 \(x_i\) 取单位根时，点值称作 DFT，插值称作 IDFT）

1. 求 \(f(x)\) 在 \(\varepsilon_{n}^{0\sim n-1}\) 的值。

   我们用分治的思路解决这个问题。但在这之前，我们假定 \(n\) 是二的幂。如果 \(n\) 不是二的幂，补若干个 \(0\) 系数使 \(n\) 是二的幂。

   \[f(x)=a_0x^0+a_1x^1+\cdots+a_{n-1}x^{n-1} \]

   定义 \(f^{(0)}(x)=a_0x^0+a_2x^1+a_4x^2+\cdots+a_{n-2}x^{(n-2)/2}\)。

   定义 \(f^{(1)}(x)=a_1x^0+a_3x^1+a_5x^3+\cdots+a_{n-1}x^{(n-2)/2}\)。

   观察到 \(f(x)=f^{(0)}(x^2)+x\cdot f^{(1)}(x^2)\)。问题转化为计算 \(f^{(0)}(x)\) 和 \(f^{(1)}(x)\) 在 \(\varepsilon_{\frac{n}{2}}^{0\sim \frac{n}{2}-1}\) 的值（求平方项的值）。

   （折半定理：\(\varepsilon_n^2=\varepsilon_{\frac{n}{2}}^{1}\)，除以 \(2\) 后再对 \(n/2\) 取模）

   于是可以递归 \(f^{(0)}(x),f^{(1)}(x)\) 求解。则 \(f(x)\) 在 \(\varepsilon_n^k\) 的点值 \(y(\varepsilon_n^k)=y^{(0)}(\varepsilon_{\frac{n}{2}}^{k\bmod \frac{n}{2}})+\varepsilon_n^ky^{(1)}(\varepsilon_{\frac{n}{2}}^{k\bmod \frac{n}{2}})\)。可以 \(O(n)\) 求得。

2. 求 \(g(x)\) 的点值。与上同理。

3. 插值还原 \(h(x)\)。也同样采取递归的思路。

   这里引入矩阵的角度理解点值。（\(\varepsilon_n^0=1\)）

   \[\begin{bmatrix} y_0\\ y_1\\ y_2\\ y_{i}\\ y_{n-1}\\ \end{bmatrix}= \begin{bmatrix} 1&1&1&\cdots&1\\ 1&\varepsilon_n^1&\varepsilon_n^2&\cdots&\varepsilon_n^{n-1}\\ 1&\varepsilon_n^2&\varepsilon_n^4&\cdots&\varepsilon_n^{2(n-1)}\\ 1&\varepsilon_n^i&\varepsilon_n^{2i}&\cdots&\varepsilon_n^{i(n-1)}\\ 1&\varepsilon_n^{n-1}&\varepsilon_n^{2(n-1)}&\cdots&\varepsilon_n^{(n-1)(n-1)}\\ \end{bmatrix} \cdot \begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1}\\ \end{bmatrix} \]

   记等号左边列向量为 \(\vec{Y}\)，矩阵为 \(V_n\)，右侧列向量为 \(\vec{A}\)。观察到 \(V_n[i][j]=\varepsilon_n^{i\cdot j}\)（编号从 \(0\) 开始）。

   \(\vec{Y}=V_n\cdot \vec{A}\)。插值就是已知 \(\vec{Y}\) 求 \(\vec{A}\)，即求 \({V_n}^{-1}\)。那 \(V_n\) 的逆 \(V_n^{-1}\) 是什么呢？

   形式也很简单：\({V_n}^{-1}[i][j]=\dfrac{1}{n}\cdot \varepsilon_n^{-i\cdot j}\)。编号也是从 \(0\) 开始。

   验证的话只要证明他俩相乘等于单位矩阵（对角线 \(1\) 其余 \(0\)）。

   观察到 \(V_n\) 和 \(V_n^{-1}\) 就是把每一个 \(\varepsilon\) 都换成了他的共轭。类似点值的递归公式 \(y(\varepsilon_n^k)=y^{(0)}(\varepsilon_{\frac{n}{2}}^{k\bmod \frac{n}{2}})+\varepsilon_n^ky^{(1)}(\varepsilon_{\frac{n}{2}}^{k\bmod \frac{n}{2}})\)，可以写出插值的递归公式 \(a(\varepsilon_n^{-k})=a^{(0)}(\varepsilon_{\frac{n}{2}}^{-k\bmod \frac{n}{2}})+\varepsilon_n^{-k}a^{(1)}(\varepsilon_{\frac{n}{2}}^{-k\bmod \frac{n}{2}})\)。

   递归的写法常数巨大。但这个过程其实可以用非递归实现。观察我们递归时哪些值被处理。

   

   （注意每次不是分首尾两半而是奇偶两半）

   发现了吗？最下面一层的二进制数如果从右往左读，刚好是 \(0\sim 2^n-1\)。

   还有一个优化：上面我们对 \(f,g\) 做了两次 DFT，对 \(h\) 做了一次 IDFT，共三次。但是如果 \(f,g\) 的系数都是实数，其实可以只做两次。

   

   摘自知乎。

至此，FFT 的算法介绍完毕。

【FFT 应用】

大整数乘积

一个 \(n\) 位的整数其实可以看作是 \(n-1\) 次多项式的 \(x\) 取 \(10\) 的结果。两个整数的乘积也可以看作是两个多项式的卷积。

合法平行四边形

题意：

考虑一个合法的平行四边形，将它补全成一个正方形。

则它的面积可以表示为 \((a+c)(b+d)-ac-bd=ad+bc=n\)。问题即求有多少个四元组 \(a,b,c,d\in\mathbb{Z_+}\) 且 \(ad+bc = n\)。

令 \(s_i\) 为 \(i\) 的因数个数。枚举 \(ad\) 和 \(bc\)，\(a,d\) 的组数就是 \(s_{ad}\)，\(b,c\) 的组数就是 \(s_{bc}\)。

所以 \(f(n)=\displaystyle\sum_{i+j=n}s_i\cdot s_j\)，即 \(f\) 是 \(s\) 的卷积。

力

把很像卷积的形式变成标准卷积：考虑拆开两半，一半标准形式，另一半对称地求。

给出 \(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_1\sim q_n\) 为 \(n\) 个排成一排的电荷。\(F_j\) 就是根据万有引力公式算出的其他电荷对 \(j\) 电荷的合力。

\(E_j=\displaystyle\sum_{i<j}\dfrac{q_j}{(i-j)^2}-\sum_{i>j}\dfrac{q_j}{(i-j)^2}\)。怎么转化为卷积的形式？

定义

\[d_k=\begin{cases}\dfrac{1}{k^2}&k>0\\0&k=0\\-\dfrac{1}{k^2}&k<0\\ \end{cases} \]

\(E_j=\sum_{i<j}q_i\cdot d_{j-i}+\sum_{i>j}q_i\cdot d_{j-i}=\sum_{i}q_i\cdot d_{j-i}\)。

但这不是标准卷积的形式，因为这里的 \(i\) 可以 \(>j\)，标准形式是不能 \(>j\) 的。这里有两种解法。

---

法一：考虑平移。

令 \(E'_{j+n}=E_j,d'_{k+n}=d_k,q_{i|i>n}=0\)。

则 \(E'_{j+n}=\sum_{i=1}^{j+n}q_i\cdot d'_{j-i+n}\)，是标准卷积形式。

---

法二：考虑对称。

令 \(U_j=\sum_{i<j}q_i\cdot d_{j-i},V_j=\sum_{i>j}q_i\cdot d_{i-j}\)。

\(U_j\) 是标准卷积形式，但是 \(V_j\) 中 \(i>j\)，不是标准形式。

考虑 \(V,q\) 进行翻转。\(V_i'=V_{n+1-i},q_i'=q_{n+1-i}\)。则 \(V_j'=\sum_{i<j}q_i\cdot d_{j-i}\)，是标准形式。

Super Rooks on Chessboard

普通的车能控制所在行列；超级车除此之外，还能控制所在的主对角线（左上-右下的对角线）。

棋盘 \(r\) 行 \(c\) 列，有 \(N\) 个超级车，超级车的位置预先给出。问：有多少个格子不被控制。

FFT + 容斥。

先反过来，求有多少个格子被至少一个控制到。

令 \(R\) 表示所有被同行的车控制的格子的集合，\(C\) 表示所有被同列的车控制的格子的集合，\(D\) 表示所有被主对角线的车控制的格子的集合。

即求 \(|R\cup C\cup D|=|R|+|C|+|D|-|R\cap C|-|R\cap D|-|C\cap D|+|R\cap C\cap D|\)。

难点在于求 \(|R\cap C\cap D|\)。

设 \(r_{1\sim n}\) 行有车，\(c_{1\sim m}\) 列有车，\(d_{1\sim l}\) 对角线有车。（对角线标号：\(y-x\)）

※ 即求有多少 \((i,j,k)\) 使得 \(c_j-r_i=d_k\)。

变形一下，\(r_i+d_k=c_j\)。考虑两个多项式 \(r(x)=\sum x^{r_i},d(x)=\sum x^{d_i}\)。

\([x^{c_j}]r(x)\cdot d(x)=\sum_{i,k}[r_i+d_k=c_j]\)。

答案就是 \(ans=\sum_i[x^{c_i}]r(x)\cdot d(x)\)。

Triple Sums

给定 \(a_1\sim a_n,a_i\le 20000,n\le 40000\)。对所有 \(S\)，问有多少三元组 \((i,j,k)\) 满足 \(i<j<k,a_i+a_j+a_k\)。

如果不限制 \(i<j<k\)。\([x^S](\sum_{i=1}^{n}x^{a_i})(\sum_{j=1}^{n}x^{a_j})(\sum_{k=1}^{n}x^{a_k})\) 是答案。

如何从上面转移到 \(\sum_{i<j<k}x^{a_i}x^{a_j}x^{a_k}\)？

\((\sum x^{a_i})^3=\sum_{i=1}^nx^{3a_i}+3\sum_{i\neq j}x^{2a_i}x^{a_j}+6\sum_{i<j<k}x^{a_i}x^{a_j}x^{a_k}\)。

要求的是 \(6\sum_{i<j<k}x^{a_i}x^{a_j}x^{a_k}\)。

等号左边用两次 FFT 可求。\(\sum_{i=1}^nx^{3a_i}\) 也很简单。真正要求的是 \(3\sum_{i\neq j}x^{2a_i}x^{a_j}\)。

考虑 \((\sum x^{2a_i})(\sum x^{a_i})=\sum_{i=1}^n x^{3a_i}+\sum_{i\neq j}x^{2a_i}x^{a_j}\)。左侧可用一次 FFT，右侧第一项很简单，于是可以求出来。

（感觉很像 \(a^2+b^2+c^2=(a+b+c)^2-2(ab+bc+ca)\) 吧？）

Point Distance

（求有哪几种距离，和有几种点对是这个距离。\(n\le 1024\)）

下面从 \(0\sim n-1\) 编号。

实际上是求 \(D_{x,y}=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}c_{i,j}\cdot c_{i+x,j+y}\)

暴力算是 \(O(n^4)\) 的。必须优化。

令 \(\phi(x,y)=x\cdot 2n+y\)。从此我们用 \(\phi\) 让一个数对应一个点对。从二维问题变成一维。这里为什么要定义成 \(\times 2n\) 而不是 \(\times n\)？因为下面 \(\phi(i,j)+\phi(x,y)\) 会用到 \(n\sim 2n\) 的位置，我们需要将它定义成 \(0\)。

\(D_{\phi(x,y)}=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}C_{\phi(i,j)}C_{\phi(i,j)+\phi(x,y)}\)。

（\(\phi(i,j)+\phi(x,y)=(i+2n+x+2n+j+y)=\phi(i+x,j+y)\)）

\(D_z=\sum_{u}C_u\cdot C_{u+z}\)，这个式子就是转化为一维的形态，对应上面那个枚举二维的式子，\(u\leftarrow \phi(i,j)\)。

（当 \(u\) 对应到 \(n-1\times n-1\) 的矩阵时，\(C_u=C_{i,j}\)；否则 \(C_u=0\)。所以有很多位置都是 \(0\)）

但是还不是卷积的形式。卷积是求和后等于左边，这里的做差后等于左边。

令 \(M=2n^2\)，令 \(C_k'=C_{M+1-K},D_K'=D_{M+1-K}\)。（翻转）

\(D'_{M+1-z}=\sum_uC_u\cdot C'_{M+1-(u+z)}\)，这下变成标准形式卷积了。

我们有 \(C\)，可以翻转求出 \(C'\)，用 FFT 把 \(C,C'\) 卷起来得到 \(D'\)，\(D'\) 翻转回去得到 \(D\)。一维的 \(D\) 转回二维后得到答案。

HDU4609

\(1\le a_1\le a_2\le\cdots\le a_n\le 1e5\)，\(n\le 1e5\)。

计数 \(i,j,k\)：\(a_i+a_j>a_k,i<j<k\)。

老样子，先忽略 \(i<j<k\) 怎么算？

对于固定的 \(k\)，答案 \(total=\sum_{w=k+1}^{2e5} [x^w](\sum_i x^{a_i})(\sum_j x^{a_j})\)。也就是多项式 \((\sum_i x^{a_i})(\sum_j x^{a_j})\) 的一个后缀和。

怎么转化？基本思路还是容斥。用 \(total\) 减去一些。

分类：

1. \(i=j\)。即 \(2a_i>a_k\)。这一种的数量也是一个后缀和 \(S\)，可以求。

2. \(i=k\neq j\)。即 \(a_i+a_j>a_i\)，这种情况肯定成立，于是这一种的 \((i,j,k)\) 的数量就是数对 \((i,j)\) 的数量 \(n\times (n-1)\)。

3. \(j=k\neq i\)。类似，这种三元组数量为 \(n\times (n-1)\)。

4. \(i,j,k\) 各不相同。

   1. \(i<j<k\) 或 \(j<i<k\)。若原题的答案记作 \(ans\)，这种的情况数是 \(2ans\)。

   2. \(k\) 不是最大值。这种情况因为 \(a\) 是升序排列的，一定满足条件。情况数就是 \(4{n\choose 3}\)。

因此 \(tot = S+ 2n\times (n-1) + 2ans + 4{n\choose 3}\).

Sum of Arrays

\(a,b\) 是长为 \(n\) 的随机生成的数列，将他们任意排列后，令 \(c_i=a_i+b_i\)，目标是让 \(c\) 中众数的出现次数最大。\(n,a_i,b_i\le 10^5\)。

先做一个简单的转化。

令 \(a[i]\) 为 \(a\) 中 \(i\) 出现的次数，\(b[i]\) 为 \(b\) 中 \(i\) 出现的次数。\(cnt[k]=\sum_{i}\min(a[i],b[k-i])\)。目标是找最大的 \(cnt[k]\)。

然后怎么转化成卷积的形式？这个 \(\min\) 我们很不喜欢。

\[\begin{aligned} cnt[k]&=\sum_{i=0}^k\sum_{x=1}^{+\infty}[a[i]\ge x][b[k-i]\ge x]\\ &=\sum_{x=1}^{+\infty}\sum_{i=0}^{k}[a[i]\ge x][b[k-i]\ge x]\\ \end{aligned} \]

令 \(\displaystyle\sum_{i=0}^k[a[i]\ge x][b[k-i]\ge x]=cnt_x[k]\)。

假定 \(x\) 固定。\(u[i]=[a[i]\ge x],v[i]=[b[i]\ge x]\)，则 \(cnt_x[k]=\displaystyle\sum_{i=0}^k u[i]\cdot v[k-i]\) 是标准形式卷积。

也就是对于固定的 \(x\)，我们可以 \(O(n\log n)\) 求出 \(cnt_x[k]\)。但是这里的 \(x\) 是枚举到 \(+\infty\) 的。

注意到题目保证 \(a,b\) 随机生成，说明 \(a[i],b[i]\) 应该不会太大。

我们可以对于 \(\displaystyle\sum_{x=1}^{20}cnt_x[k]\) 用 FFT 求，复杂度 \(O(20n\log n)\)；对于 \(\displaystyle\sum_{x=21}^{+\infty}\) 的，可以直接暴力找出所有 \(>21\) 的 \(a[i],b[i]\)。虽然复杂度玄学，但因为随机数据，跑得很快。

【NTT 算法】

NTT，数论变换。它解决的问题是：两个多项式 \(A,B\) 之积，结果每一项的系数都对质数 \(p\) 取模。

如果用 FFT 最后再取模有个问题：FFT 算出来的大多数都是浮点数，精度不足。

众所周知，FFT 是先在 \(n\) 次单位根的位置点值，然后再插值。（这里先把 \(A,B\) 扩充成 \(2\) 的幂）

阶：\(a\) 模 \(b\) 的阶记作 \(\delta_b(a)\)，\(\delta_b(a)\) 是所有满足 \(a^x\equiv 1\pmod b\) 的正整数 \(x\) 中最小的那个。

原根：若 \(\delta_b(a)=\phi(b)\)，称 \(a\) 是模 \(b\) 的原根。