FFT和NTT学习笔记
算法导论
http://picks.logdown.com/posts/177631-fast-fourier-transform
https://blog.csdn.net/qq_38944163/article/details/81835205
https://www.cnblogs.com/RabbitHu/p/FFT.html
概述
目的
以(O(nlg_n))的时间复杂度计算多项式乘法
多项式的表达
- 系数表达: ({a_0, a_1, ..., a_{n-1}})
- 点值表达: ({(x_0,y_0), (x_1,y_1), ..., (x_{n-1},y_{n-1})})
- 插值多项式唯一性: (x)各不相同得到唯一的多项式
策略
关于单位复数根
复平面
复数平面（complex plane）是用水平的实轴与垂直的虚轴建立起来的复数的几何表示。
一个复数的实部用沿着 x-轴 的位移表示，虚部用沿着 y-轴 的位移表示。
复数乘法： 模长相乘， 幅角相加（幅角：向量与x轴正向夹角）
单位复数根
感性理解:
将复平面上的单位圆从 x 轴起分成 n 等分， 从 x 轴起标号 0..n-1
那么由 “模长相乘， 幅角相加” 可得， 标号为 1 的那个复数的 (k) 次方， 就是标号为 (k) 的
所以将这些数记为 (omega_n^0, omega_n^1, dots, omega_n^{n-1})
理性证明:
- (n)次单位复数根 是满足 (omega^n=1) 的复数 (omega), 恰好有 (n) 个
由 (e^{i heta}=cos( heta)+isin( heta))
得 (e^{2pi i}=cos(2pi)+isin(2pi)=1)
- 所以对于 (k=0,1,...n-1) , 这些根是 (e^{2pi i k/n})
- 有复数乘法可得这 (n) 个单位复数根均匀分布在以复平面原点为圆心的单位半径的圆周上
- (omega_n=e^{2pi i/n}) 称为 主(n)次单位根，其他单位根都是 (omega_n) 的幂次，并有 (omega_n^{j+k}=omega_n^{(j+k) mod n})
单位复数根基本性质
(消去引理)
- 对任何整数 (以及ngeq0,kgeq0,以及dgeq0) , (omega_{dn}^{dk}=omega_n^k)
- 感性理解： 原来分成 (n) 等分， 取第 (k) 个， 和分成 (dn) 等分， 取第 (dk) 个当然一样
- 证明：根据指数形式定义式
(折半引理)
感性理解：
(omega_n^{k + n/2} = -omega_n^k)
复平面单位圆上关于原点对称
理性证明：
- 如果 (n>0) 为偶数，那么 (n) 个 (n) 次单位复数根的平方的集合就是 (n/2) 个 (n/2) 次单位复数根的集合
- 证明：((omega_n^{k+n/2})^2=omega_n^{2k+n}=omega_n^{2k}omega_n^n=(omega_n^k)^2=omega_{n/2}^k)
- 也可以根据 (omega_n^{n/2}=-1) 得到 (omega_n^{k+n/2}=-omega_n^k)
- 由此可见他可以递归让子问题的规模缩小一半
(求和引理)
感性理解：
(n) 是 (2) 的次幂
(k=0) 显然
(k
eq 0) 时
首先一个偶数 (n) 的等分点的和是 (0) (对称点), 那么对于 (n=2^k) 显然适用
假如 (k) 是 (n) 的因数, 那么 (k) 也是 (2) 的因数, 它会取一个等分的循环取若干次, 和是 (0)
否则互质, 那么它会取遍这 (n) 个点, 和也是 (0)
- 证明：等比数列求和 (sum_{j=0}^{n-1}(omega_n^k)^j=frac{1*(1-omega_n^{kn})}{1-omega_n^k}=0)
- 因为要求 (k) 不能被 (n) 整除，保证分母不为 (0)
概述
系数转点值: 求出 (n) 个单位根的点值 ((y_0, y_1, dots, y_{n-1}))
点值转系数: 将这 (n) 个点值作为一个新的多项式 (B(x)) 的系数, 用单位根的倒数求一次点值 ((z_0, z_1, dots, z_{n-1}))
展开可得
根据之前的求和引理可得: (z_k = n cdot a_k)
于是系数 (a_k = z_k / n)
也就是说, 这两个过程都可以通过求点值来完成,
所以要完成两个多项式的相乘, 只需要先求一遍他们的点值, 点值相乘, 在转成系数就行了
DFT
点值向量(y=(y_0, y_1, ..., y_{n-1}))就是系数向量(a=(a_0, a_1, ..., a_{n-1}))的离散傅里叶变换(DFT), 也记为(y=DFT_n(a))
FFT
快速傅里叶变换(FFT)，利用单位复数根的特殊性质，在(O(nlg_n))的时间内计算出(DFT_n(a))
分治策略，采用(A(x))中偶数下标和奇数下标的系数，分别定义两个次数界为(n/2)的多项式
- (A_0(x) = a0 + a2*x + ... + a_{n-2}*x^{n/2-1})
- (A_1(x) = a1 + a3*x + ... + a_{n-1}*x^{n/2-1})
(A(x) = A_0(x^2) + x cdot A_1(x^2))
带入单位根可以发现
这两个关于远点对称的单位根的点值只是相差一个符号, 所以计算 (A_0(omega_{frac{n}{2}}^{k}), omega_n^kA_1(omega_{frac{n}{2}}^{k})) 就可得到两个单位根的取值
因为应用了(omega_n^k)的正负数形式，所以称其为旋转因子
可得到以下伪代码，时间复杂度(O(nlg_n))
RECURSIVE_FFT(a) // 计算DFTn(a)
n = a.length
if n == 1
return a // a * (x^0)
a[0] = (a(0), a(2), ..., a(n-2))
a[1] = (a(1), a(3), ..., a(n-1))
y[0] = RECURSIVE_FFT(a[0])
y[1] = RECURSIVE_FFT(a[1])
for k = 0 to n/2-1
y[k] = y[0][k] + omega[n][k] * y[1][k]
y[k+n/2] = y[0][k] - omega[n][k] * y[1][k]
return y;
IDFT
通过点值得到系数
IFFT
在概述中说明了, 只需要将点值作为系数的到新的多项式, 再带入单位根的倒数, 得到这个多项式的点值
然后点值 (/n) 即可得到系数
高效FFT实现
画出系数的递归树可得叶子的标号的规律:
蝴蝶操作中
可以观察到(y_k^{[1]})的位置其实就是(y_{k+n/2})的位置，相当于两个位置上的(y)值进行蝴蝶操作并只对这两个位置进行修改
观察下图递归FFT输入的向量树
其叶子有这样的规律：
记(rev(x))为(x)的二进制表示的倒串所形成的(lg_n)位的数,则(rev(x))位上为(a_x)
因为每次的叶子是一样的( (n) 相同), 所以可以预处理, 并且可以 (O(n)) 递推 rev[]
注意 预处理叶子位置不是扫到 (len / 2)
不妨从下往上模拟递归的过程
方便起见(省去二维的数组), 将每一层中各多项式的各点值标号, 将递归的写成这样
{
FFT(a, n / 2);
FFT(a + n / 2, n / 2);
for (int i = 0; i < n / 2; ++ i) // 枚举 omega_n^i
{
buf[i] = a[i] + a[i + m];
buf[i + m] = a[i] - omega_{n / 2}^{-i} * a[i + m];
}
}
转成非递归, 现将原系数变成叶子的状态
for (int l = 2; l <= len; l <<= 1)
{
int ll = l >> 1;
for (int i = 0; i < len; i += l)
{
for (int j = 0; j < l / 2; ++ j)
{
CPD x = ome[len / l * j];
if (opt) x = inv[len / l * j];
buf[i + j] = a[i + j] + x * a[i + j + ll];
buf[i + j + ll] = a[i + j] - x * a[i + j + ll];
}
}
for (int i = 0; i < len; ++ i) a[i] = buf[i];
}
注意这里的单位根可以利用引理对应到 (n) 次单位根上
卡常提示:
- 然后注意不要预处理
ome[], inv[], 循环的时候直接乘上第一个单位根即可 - 手写复数, 并且不要用构造函数(常数小了一半)
- 2 * PI / l 可以改为 PI / mid
// P1919 A*B(FFT)
//#pragma GCC optimize(2)
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
struct Complex
{
double x, y;
Complex operator + (const Complex &b) const { return (Complex){x + b.x, y + b.y}; }
Complex operator - (const Complex &b) const { return (Complex){x - b.x, y - b.y}; }
Complex operator * (const Complex &b) const { return (Complex){x * b.x - y * b.y, x * b.y + y * b.x}; }
} ;
const int MAXN = (1 << 21) + 3;
char ch[MAXN];
int lena, lenb, n, dgt;
Complex a[MAXN], b[MAXN];
void read(Complex * a, int & len)
{
scanf("%s", ch + 1); len = strlen(ch + 1);
for (int i = 0; i < len; ++ i) a[i].x = ch[len - i] - '0';
}
int rev[MAXN];
void init(int n, int dgt)
{
for (int i = 0; i < n; ++ i)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (dgt - 1));
}
const double PI = acos(-1.0);
void FFT(Complex * a, int len, int opt)
{
for (int i = 0; i < len; ++ i)
if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int mid = 1; mid < len; mid <<= 1)
{
Complex omen = (Complex){cos(PI / mid), opt * sin(PI / mid)} ;
for (int i = 0; i < len; i += (mid << 1))
{
Complex ome = (Complex){1, 0} ;
for (int j = 0; j < mid; ++ j, ome = ome * omen)
{
Complex t = ome * a[i + j + mid];
a[i + j + mid] = a[i + j] - t;
a[i + j] = a[i + j] + t;
}
}
}
}
int ret[MAXN];
/*
20191212
0724~0757~0839
a * b FFT
*/
int main()
{
read(a, lena); read(b, lenb);
n = 1; dgt = 0;
while (n < lena + lenb) n <<= 1, ++ dgt;
init(n, dgt);
FFT(a, n, 1);
FFT(b, n, 1);
for (int i = 0; i < n; ++ i) a[i] = a[i] * b[i];
FFT(a, n, -1);
for (int i = 0; i < n; ++ i)
{
ret[i] += (int)(a[i].x / n + 0.5);
ret[i + 1] += ret[i] / 10;
ret[i] %= 10;
}
while (!ret[n - 1] && n > 0) -- n;
for (int i = n - 1; i >= 0; -- i) printf("%d", ret[i]);
// printf("#%d
", cnt);
return 0;
}
/*
76543210
76543210
6543
21
*/
总结
FFT 利用了单位根的以下性质
from http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform#i-17
首先来看 FFT 中能在 O(nlogn) 时间内变换用到了单位根 ω 的什么性质 ...
自己解释一遍:
- 单位根 (omega_n^0, dots, omega_n^{n-1}) 是单位圆的等分点(不同), 所以有了 "消去引理"
- (omega_n^{frac{n}{2}+k}=-omega_n^k, omega_n^2=omega_{frac{n}{2}})
- 根据奇偶系数拆分的表达式, 这一层的 (m) 个多项式都要带入 (n) 个值, 可以通过下一层 (2m) 个多项式带入 (n/2) 个值求得, 由此可知 (log) 层
- (sum_{j=0}^{n-1}(omega_n^k)^j = egin{cases} n, & k = 0 \ 0, & k eq 0 \ end{cases}) 使得点值转系数的逆变换可以同样方式完成
其中第一点是最基本的, 他推出了后面的性质
NTT
在取一些模数的意义下, 原根具有上述性质
区别在于原根的手动计算不同
void NTT(LL * a, int len, int sgn)
{
for (int i = 0; i < len; ++ i)
if (i < rev[i]) swap(a[i], a[rev[i]]);
for (int mid = 1; mid < len; mid <<= 1)
{
for (int i = 0; i < len; i += (mid << 1))
{
for (int j = 0; j < mid; ++ j)
{
LL t = (sgn == 1 ? g[len / mid / 2 * j] : ig[len / mid / 2 * j]) * a[i + j + mid] % MOD;
a[i + j + mid] = (a[i + j] + MOD - t) % MOD;
a[i + j] = (a[i + j] + t) % MOD;
}
}
}
}
invn = inv(n);
g[0] = ig[0] = 1;
g[1] = fpow(3, (MOD - 1) / n); ig[1] = inv(g[1]);
for (int i = 2; i < n; ++ i)
g[i] = g[i - 1] * g[1] % MOD, ig[i] = ig[i - 1] * ig[1] % MOD;