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)
-
求 \(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)\) 求得。
-
求 \(g(x)\) 的点值。与上同理。
-
插值还原 \(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\),定义
对 \(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}\)。怎么转化为卷积的形式?
定义
\(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\) 减去一些。
分类:
-
\(i=j\)。即 \(2a_i>a_k\)。这一种的数量也是一个后缀和 \(S\),可以求。
-
\(i=k\neq j\)。即 \(a_i+a_j>a_i\),这种情况肯定成立,于是这一种的 \((i,j,k)\) 的数量就是数对 \((i,j)\) 的数量 \(n\times (n-1)\)。
-
\(j=k\neq i\)。类似,这种三元组数量为 \(n\times (n-1)\)。
-
\(i,j,k\) 各不相同。
-
\(i<j<k\) 或 \(j<i<k\)。若原题的答案记作 \(ans\),这种的情况数是 \(2ans\)。
-
\(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\) 我们很不喜欢。
令 \(\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\) 的原根。