FFT&NTT总结
FFT&NTT总结
一些概念
\(DFT:\)离散傅里叶变换\(\rightarrow O(n^2)\)计算多项式卷积
\(FFT:\)快速傅里叶变换\(\rightarrow O(nlogn)\)计算多项式卷积
\(NTT:\)快速数论变换\(\rightarrow\)对\(FFT\)的常数优化
\(MTT:\)\(NTT\)的一些拓展
FFT
多项式&卷积
设\(A(x)\)表示一个\(n-1\)次多项式
则\(A(x)=\sum_{i=0}^{n-1}a_ix^i\)
而卷积就是两个多项式相乘
如果我们像平常那样暴力乘起来,复杂度是\(O(n^2)\)的
点值表示法
将\(n\)个点代进\(n-1\)次多项式\(A(x)\)
则可以确定\(n\)对\((x,y)\)
而我们有\(n\)个点也可以确定一个\(n-1\)次多项式
为什么?很(bu)显(hui)然(zheng)啊。
我们后面的\(FFT\)的优化就是基于这个来的
复数
定义
我们把形如\(z=a+bi\)(\(a,b\)均为实数)的数称为复数,其中\(a\)称为实部,\(b\)称为虚部,\(i\)称为虚数单位。(摘自百度百科)其中\(i^2=-1\)。
而在复平面中,\(x\)轴代表实数,\(y\)轴代表虚数(除原点),从原点\((0,0)\)到\((a,b)\)代表复数\(a+bi\)
模长:\((0,0)\)到\((a,b)\)的距离,即\(\sqrt {a^2+b^2}\)
幅角:以逆时针为正方向,\(x\)轴到已知向量的转角的有向角
运算法则
加减法:
和向量一样,即
\((a,b)+(c,d)=(a+b,c+d)\)
\((a,b)-(c,d)=(a-b,c-d)\)
乘法:
几何意义:复数相乘,模长相乘,幅角相加
代数定义:
单位根
(下文默认\(n\)为\(2\)的整数次幂)
在复平面上,以原点为圆心,\(1\)为半径的圆叫做单位圆。
以原点为起点,圆的\(n\)等分点为终点,作\(n\)个向量,设幅角为正且最小的复数向量为\(\omega _n\),称为\(n\)次单位根。
\(n\)个向量为\(\omega_n^1,\omega_n^2,\omega_n^3...\omega_n^{n-1},\omega_n^n\)(\(\omega_n^n=\omega_n^0=1\))
如何计算他们的值呢,
可以用欧拉公式:
单位根的幅角为周角的\(\frac 1n\)
代数中,若\(z^n=1\),我们把\(z\)称为\(n\)次单位根
性质
1、\(\omega_n^k=cos\;(k*\frac{2\pi}{n})+i*sin\;(k*\frac{2\pi}{n})\)
2、\(\omega_n^k=\omega_{2n}^{2k}\)
3、\(\omega_n^{k+\frac n2}=-\omega_n^k\)
4、\(\omega_n^0=\omega_n^n=1\)
快速傅里叶变换
我们前面提过,一个\(n-1\)次多项式可以用\(n\)个点唯一确定,
我们可以把\(0\)~\(n-1\)次单位根依次带入
但仍然是\(O(n^2)\)啊,因为单位根有很多优秀的性质
所以我们来推一波公式
有
按照下表奇偶性分类
设
则
将\(\omega_n^k(k<\frac n2)\)代入\(:A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\)
将\(\omega_n^{k+\frac n2}\)代入:\(A(\omega_n^{k+\frac n2})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\)
发现只有一个符号不一样
于是求第一个式子时,我们可以\(O(1)\)求第二个式子
我们就将这个问题缩小了一半
递归搞下去,就可以\(O(nlogn)\)了
快速傅里叶逆变换
真的不想写了2333
跟上面其实差不多,直接看代码吧。。。
下面代码中
FFT(a, -1);
for (int i = 0; i <= M; i++) printf("%d ", (int)(a[i].x / N + 0.5));
是快速傅里叶逆变换
代码
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
namespace IO {
const int BUFSIZE = 1 << 20;
char ibuf[BUFSIZE], *is = ibuf, *it = ibuf;
inline char gc() {
if (is == it) it = (is = ibuf) + fread(ibuf, 1, BUFSIZE, stdin);
return *is++;
}
}
inline int gi() {
register int data = 0, w = 1;
register char ch = 0;
while (!isdigit(ch) && ch != '-') ch = IO::gc();
if (ch == '-') w = -1, ch = IO::gc();
while (isdigit(ch)) data = 10 * data + ch - '0', ch = IO::gc();
return w * data;
}
const double PI = acos(-1.0);
const int MAX_N = 3e6 + 5;
struct Complex { double x, y; } a[MAX_N], b[MAX_N];
Complex operator + (const Complex &a, const Complex &b) { return (Complex){a.x + b.x, a.y + b.y}; }
Complex operator - (const Complex &a, const Complex &b) { return (Complex){a.x - b.x, a.y - b.y}; }
Complex operator * (const Complex &a, const Complex &b) { return (Complex){a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x}; }
int N, M, P, r[MAX_N];
void FFT(Complex *p, int op) {
for (int i = 0; i < N; i++) if (i < r[i]) swap(p[i], p[r[i]]);
for (int i = 1; i < N; i <<= 1) {
Complex rot = (Complex){cos(PI / i), op * sin(PI / i)};
for (int j = 0; j < N; j += (i << 1)) {
Complex w = (Complex){1, 0};
for (int k = 0; k < i; ++k, w = w * rot) {
Complex x = p[j + k], y = w * p[j + k + i];
p[j + k] = x + y, p[j + k + i] = x - y;
}
}
}
}
int main () {
N = gi(), M = gi();
for (int i = 0; i <= N; i++) a[i].x = gi();
for (int i = 0; i <= M; i++) b[i].x = gi();
for (M += N, N = 1; N <= M; N <<= 1, ++P) ;
for (int i = 0; i < N; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (P - 1));
FFT(a, 1), FFT(b, 1);
for (int i = 0; i < N; i++) a[i] = a[i] * b[i];
FFT(a, -1);
for (int i = 0; i <= M; i++) printf("%d ", (int)(a[i].x / N + 0.5));
return 0;
}
NTT
其实和\(FFT\)差不多啦,
就是把单位根换为原根就行了
代码
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;
inline int gi() {
register int data = 0, w = 1;
register char ch = 0;
while (!isdigit(ch) && ch != '-') ch = getchar();
if (ch == '-') w = -1, ch = getchar();
while (isdigit(ch)) data = 10 * data + ch - '0', ch = getchar();
return w * data;
}
const int MAX_N = 3e6 + 5, Mod = 998244353, G = 3, iG = 332748118;
int fpow(int x, int y) {
int res = 1;
while (y) {
if (y & 1) res = 1ll * res * x % Mod;
x = 1ll * x * x % Mod;
y >>= 1;
}
return res;
}
int Limit = 1, r[MAX_N];
void NTT(int *p, int op) {
for (int i = 0; i < Limit; i++) if (i < r[i]) swap(p[i], p[r[i]]);
for (int i = 1; i < Limit; i <<= 1) {
int rot = fpow(op == 1 ? G : iG, (Mod - 1) / (i << 1));
for (int j = 0, pls = (i << 1); j < Limit; j += pls) {
int w = 1;
for (int k = 0; k < i; k++, w = 1ll * w * rot % Mod) {
int x = p[j + k], y = 1ll * w * p[i + k + j] % Mod;
p[j + k] = (x + y) % Mod, p[i + j + k] = (x - y + Mod) % Mod;
}
}
}
}
int N, M, a[MAX_N], b[MAX_N];
int main () {
N = gi(), M = gi();
for (int i = 0; i <= N; i++) a[i] = (gi() + Mod) % Mod;
for (int i = 0; i <= M; i++) b[i] = (gi() + Mod) % Mod;
int L = 0;
while (Limit <= N + M) Limit <<= 1, ++L;
for (int i = 0; i < Limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (L - 1));
NTT(a, 1), NTT(b, 1);
for (int i = 0; i < Limit; i++) a[i] = 1ll * a[i] * b[i] % Mod;
NTT(a, -1);
int inv = fpow(Limit, Mod - 2);
for (int i = 0; i <= N + M; i++) printf("%lld ", (1ll * a[i] * inv) % Mod);
return 0;
}