模板及讲解
简介
$DFT$：离散傅里叶变换
$IDFT$：离散傅里叶逆变换
$FFT$：快速傅里叶变换，计算多项式乘法 (卷积)
$FNTT/NTT$：快速傅里叶变换的优化版，优化常数及误差
FFT
多项式
系数表示法
设$A(x)$表示一个$n−1$次多项式
则$A(x)=\sum\limits_{i=0}^n a_ix^i$
利用这种方法计算多项式乘法复杂度为$O(n^2)$
（第一个多项式中每个系数都需要与第二个多项式的每个系数相乘）
点值表示法
将$n$互不相同的$x$带入多项式，会得到$n$个不同的取值$y$
则该多项式被这$n$个点$(x_1,y_1),(x_2,y_2),\dots,(x_n,y_n)$唯一确定
其中$y_i=\sum\limits_{j=0}^{n-1} a_jx_i^j$
利用这种方法计算多项式乘法复杂度为$O(n^2)$
点值表示法用处
我们可以把系数表示转化成点值表示。点值表示下的多项式就是选相同的$x_i$，把对应的$y_i$相乘。
两个长度为$n$的多项式相乘得到的是长度为$2n−1$的多项式，需要$2n−1$个点值才能唯一表示，因此一开始两个多项式也要选$2n−1$个点表示。
随意选$O(n)$个点，计算$y$是$O(n)$的，总时间还是$O(n^2)$的，所以我们要介绍后面的利器，单位根
复数
定义
设$a,b$为实数，$i^2=−1$，形如$a+bi$的数叫复数，其中$i$被称为虚数单位
在复平面中，$x$代表实数，$y$轴（除原点外的点）代表虚数，从原点$(0,0)$到$(a,b)$的向量表示复数$a+bi$
模长：从原点$(0,0)$到点$(a,b)$的距离，即$\sqrt{a^2+b^2}$
幅角：假设以逆时针为正方向，从$x$轴正半轴到已知向量的转角的有向角叫做幅角
运算法则
加法：$(a,b)+(x,y)=(a+x,b+y)$
乘法：$(a,b) \times (x,y)=(ac-bd, bc+ad)$ (可从定义式中推出)
单位根
定义
$n$次单位根（记为$\omega$）
它是$n$个复数的集合，每一个的$n$次方都等于$1$。其中的一个是$e^{\frac {2\pi i} n}$，记为$\omega_n$。
欧拉公式
$$
e^{ix}=\cos x+(\sin x)i
$$
$x$就是这个复数向量的旋转角。
易得
$$
\omega_{n}^{k}=\cos k \frac{2\pi}{n}+i\sin k\frac{2\pi}{n}
$$
几何意义
在复平面上，以原点为圆心，$1$为半径作圆，所得的圆叫单位圆。以圆点为起点，圆的$n$等分点为终点，做$n$个向量，设幅角为正且最小的向量对应的复数为$\omega_n$，称为$n$次单位根。
下面是算导的图片，这$n$个复数向量在单位圆上呈放射状。
性质
1、消去引理 (从定义式可推，可从图中观察)：
$$
\omega _{2n}^{2k}=\omega _{n}^{k}
$$
2、(从定义式可推，可从图中观察)
$$
\omega _{n}^{k+\frac{n}{2}}=-\omega _{n}^{k}
$$
3、(从定义式可推，可从图中观察)
$$
\omega_n^0=\omega_n^n=1
$$
4、(可从图中观察)
$$
\omega^{\frac n 2}_n=-1
$$
$DFT$和$IDFT$
$DFT$
我们把$\omega ^k_n(k \in [n−1])$分别带入多项式求点值
假设$n$为偶数，那么我们可以把它的奇偶项分开，用两个多项式表示它
$$
A^{[0]}(x)=a_0+a_2x+a_4x^2+…+a_{n−2}x^{\frac n2−1} \\
A^{[1]}(x)=a_1+a_3x+a_5x^2+…+a_{n−1}x^{\frac n2−1}
$$
则 $A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$
带入单位根$((k < \frac n 2))$
$$
A(\omega^k_n)=A^{[0]}(\omega^{2k}_n)+\omega^k_nA^{[1]}(\omega^{2k}_n) \\
=A^{[0]}(\omega^k_{\frac n 2})+\omega^k_nA^{[1]}(\omega^k_{\frac n 2})(k < \frac n 2)
$$
$((k \geq \frac n 2))$时
$$
A(\omega^{\frac n 2+k}_n)=A^{[0]}(\omega^{n+2k}_n)+\omega^{\frac n 2+k}_nA^{[1]}(\omega^{n+2k}_n)\\
=A^{[0]}(\omega^{2k}_n)+\omega^{\frac n 2}_n\omega^k_nA^{[1]}(\omega^{2k}_n) \\
=A^{[0]}(\omega^k_{\frac n 2})-\omega^k_nA^{[1]}(\omega^k_{\frac n 2})
$$
可以由算法写出$DFT$的复杂度$T(n)=2T(\frac n 2)+O(n)=O(n\log n)$
$IDFT$
还原成系数表示的过程叫做$IDFT$。
只需要记住，$IDFT$的$\omega_n$是$e^{-\frac{2\pi i}n}$，最后的结果除以$n$，其它过程同$DFT$，可以写在一个函数里。
蝴蝶操作和$Rader$排序
在操作过程中，取出了a0[k]和a1[k]的值，通过计算又求出了a[k]和a[k+n]的值。我们把这样的一次运算叫做「蝴蝶操作」。
以$n=8$为例，看看递归过程的结构
我们完全可以从底向上递推求解。
剩下的问题就是把初始的数组变成最后一层的样子了。
观察最后序列的编号的二进制表示000,100,010,110,001,101,011,111与原来000,001,010,011,100,101,110,111相比，每个位置上的二进制位都反过来了。这样的变化叫做$Rader$排序。
我们可以$O(n)$将$Rader$排序的映射关系求出。设$iRader$排序后的数为$r_i$，我们可以通过$r_{\frac i2}$递推求出
这个递推式相当于偶数位$i$的值就是$\frac i2$值的一半，而奇数位则在前一个偶数位上$\operatorname{or} 2^{l-1}$，其中$l$是$2^l$超过多项式最终长度的最小$l$
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<map>
#include<set>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db long double
#define fir first
#define sec second
#define mp make_pair
using namespace std;
namespace flyinthesky {
const int MAXN = 4200000 + 5;
const double PI = acos(-1);
struct C { // 手写 complex 类
double r, i;
C() {r = i = 0;}
C(double x, double y){r = x, i = y;}
C operator + (C &x) {return C(r + x.r, i + x.i);}
C operator - (C &x) {return C(r - x.r, i - x.i);}
C operator * (C &x) {return C(r * x.r - i * x.i, r * x.i + i * x.r);}
void operator += (C &x) {r += x.r, i += x.i;}
void operator *= (C &x) {
double t = r;
r = r * x.r - i * x.i, i = t * x.i + i * x.r;
}
}f[MAXN], g[MAXN];
int n, m, r[MAXN];
void FFT(C *a, int op) { // op = 1 DFT; op = -1 IDFT;
for (int i = 0; i < n; ++i) // 根据映射关系交换元素
if (i < r[i]) swap(a[i], a[r[i]]);
for (int i = 1; i < n; i *= 2) { // 第几层
C Wn = C(cos(PI / i), sin(PI / i) * op);
for (int j = 0; j < n; j += i * 2) { // 一层中的子问题 [j, j + i * 2 - 1]
C w = C(1, 0), *a0 = a + j, *a1 = a0 + i;
// a0 = [j, j + i - 1], a1 = [j + i, j + i * 2 - 1]
for (int k = 0; k < i; ++k) { // 处理子问题
C t = *a1 * w;
*a1 = *a0 - t, *a0 += t; // 蝴蝶操作
// a1 = a0 - a1w
// a0 = a0 + a1w
w *= Wn, ++a0, ++a1;
}
}
}
}
void clean() {
}
int solve() {
clean();
scanf("%d%d", &n, &m);
for (int x, i = 0; i <= n; ++i) scanf("%d", &x), f[i].r = (db)x;
for (int x, i = 0; i <= m; ++i) scanf("%d", &x), g[i].r = (db)x;
int l = 0;
for (m += n, n = 1; n <= m; n <<= 1, ++l);
for (int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); //递推求r
FFT(f, 1), FFT(g, 1);
for (int i = 0; i < n; ++i) f[i] *= g[i]; // 对应 y 相乘
FFT(f, -1);
for (int i = 0; i <= m; ++i) printf("%.0f ",fabs(f[i].r) / n); // IDFT 要除以 n
return 0;
}
}
int main() {
flyinthesky::solve();
return 0;
}
/*
1 6
5 6
1 2 3 4 5 6 7
*/
写时注意
1、fabs(f[i].r)
2、最好四舍五入
NTT
FFT运用复数会有误差，并且不能取模，所以就有了NTT
用原根代替单位根即可。
原根
如果$G^n=1(\mod p)$中$n$的最小正整数解为欧拉函数$\varphi(p)$，则称$G$为$p$的一个原根。
因为$G^{p−1}=\omega^n_n=1$，所以有$\omega_n=g^{\frac{p−1}n}$，并且原根也满足原来FFT用的单位根的性质
常用NTT模数：$998244353(998244352=2^{23}×7×17)$，$3$是它的一个原根
代码
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cmath>
#include<map>
#include<set>
#define ms(i, j) memset(i, j, sizeof i)
#define LL long long
#define db long double
#define fir first
#define sec second
#define mp make_pair
using namespace std;
namespace flyinthesky {
const LL MAXN = 4200000 + 5, MO = 998244353, G = 3;
LL f[MAXN], g[MAXN];
LL n, m, r[MAXN];
LL ksm(LL a, LL b) {
LL bs = a, ans = 1;
while (b) {
if (b & 1) ans = ans * bs % MO;
bs = bs * bs % MO;
b >>= 1;
}
return ans;
}
void NTT(LL *a, LL op) { // op = 1 DFT; op = -1 IDFT;
for (LL i = 0; i < n; ++i) // 根据映射关系交换元素
if (i < r[i]) swap(a[i], a[r[i]]);
for (LL i = 1; i < n; i <<= 1) { // 第几层
LL Wn = ksm(G, (MO - 1) / (i << 1));
if (op == -1) Wn = ksm(Wn, MO - 2);
for (LL j = 0; j < n; j += (i << 1)) { // 一层中的子问题 [j, j + (i << 1) - 1]
LL w = 1, *a0 = a + j, *a1 = a0 + i;
// a0 = [j, j + i - 1], a1 = [j + i, j + (i << 1) - 1]
for (LL k = 0; k < i; ++k) { // 处理子问题
LL t = *a1 * w % MO;
*a1 = (*a0 - t + MO) % MO, *a0 = (*a0 + t) % MO; // 蝴蝶操作
// a1 = a0 - a1w
// a0 = a0 + a1w
w = w * Wn % MO, ++a0, ++a1;
}
}
}
}
void clean() {
}
int solve() {
clean();
scanf("%lld%lld", &n, &m);
for (LL i = 0; i <= n; ++i) scanf("%lld", &f[i]);
for (LL i = 0; i <= m; ++i) scanf("%lld", &g[i]);
LL l = 0;
for (m += n, n = 1; n <= m; n <<= 1, ++l);
for (LL i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); //递推求r
NTT(f, 1), NTT(g, 1);
for (LL i = 0; i < n; ++i) f[i] = f[i] * g[i] % MO; // 对应 y 相乘
NTT(f, -1);
LL inv = ksm(n, MO - 2);
for (LL i = 0; i <= m; ++i) printf("%lld ", f[i] * inv % MO); // IDFT 要除以 n
return 0;
}
}
int main() {
flyinthesky::solve();
return 0;
}
/*
1 6
5 6
1 2 3 4 5 6 7
*/
卷积
定义
卷积就是形如
$$
C_n=\sum\limits_{i=0}^{n}a_ib_{n-i}
$$
卷积的两项底数相加要为一个常数。
乘法竖式就是一个卷积，所以我们可以用 FFT 求大整数乘法
上述多项式乘法也是卷积
