FFT（1）
FFT 是一种实现两个 \(n\) 项多项式相乘的 \(O(n\log n)\) 高效算法。它可以用来快速求解大数乘法，还在许多其他方面有广泛应用。
FFT 的介绍文章大多数逼格较高，让没有一定高中数学基础的同学难以有耐心阅读下去、接受。当然，胡小兔：小学生都能看懂的FFT 是非常便于理解的。本文面向对象，致力于降低接受难度，但不一定能达到目标。
定义
FFT 是实现两个 \(n\) 项多项式 \(A(x)=a_0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1},B(x)=b_0+b_1x^1+b_2x^2+...+b_{n-1}x^{n-1}\) 相乘的 \(O(n\log n)\) 快速算法。
下文规定 \(n\) 为 2 的整数幂（并且 FFT 也仅能处理项数为 2 的整数幂的多项式）。但通常 \(n\) 不会是 2 的整数幂，这时在多项式后补 0，即 \(A(x)=a_0+a_1x^1+...+0\cdot x^{n-2}+0\cdot x^{n-1}\)。
多项式的系数表示和点值表示
通常我们看到的多项式写成 \(A(x)=a_0+a_1x^1+a_2x^2+...+a_{n-1}x^{n-1}\) 就是系数表示。
引理 一个 \(n-1\) 次多项式 \(A(x)\) 可以被 \(n\) 个点 \((x_0,A(x_0)),(x_1,A(x_1)),...,(x_{n-1},A(x_{n-1}))\) 唯一确定。（\(x_0,x_1,...,x_{n-1}\) 互不相等。）例如，两点确定一条直线，三点确定一条二次函数。这是因为一个 \(n\) 个方程的方程组可以解出 \(a_0\sim a_{n-1}\) 共 \(n\) 个参数（未知数）。
因此
用自变量和函数值的坐标系内点表达一个多项式的方法就是点值表示。
点值表示有什么用？
第一点，点值表示下，两个多项式相乘是 \(O(n)\) 的：
假如我们可以实现系数表示转为点值表示和点值表示转为系数表示，就可以完成我们的问题。
目前来看，系数表示转为点值表示我们需要 \(O(n^2)\)，点值表示转为系数表示需要 \(O(n^3)\) 高斯消元。FFT 有更高效方法。
单位复根
这个名字听起来十分的可怕，但其实理解其思想也没有那么难。
在了解单位复根之前，我们首先要了解复数和向量。
什么是向量？
向量可以理解为坐标平面内一条具有方向和长度的线段。（并且它表达这两种属性，没有所谓“位置”。）它的长度叫做「模长」。
同时，我们发现如果把线段的箭尾与原点重合，便可以用箭头的坐标来唯一确定一个向量。以向量作为终边的角称作「极角」。
向量乘法：模长相乘，极角相加。
什么是复数？
我们都听说过 \(\mathbf{i}=\sqrt{-1}\)。复数就是由所有的 \(a+b\mathbf{i}(a,b\in \R)\) 构成的数集。复数分为实数（\(b=0\)）和虚数（\(b\ne 0\)），虚数又有纯虚数（\(a=0\) 且 \(b\ne 0\)）。
复数 \(z=a+b\mathbf{i}\) 可以看成坐标平面内的点 \((a,b)\)。下面称这个坐标平面为￥。
根据复数的性质，数学家傅里叶找到了一种方法，使得带入一系列特殊的 \(x_0,x_1,...,x_{n-1}\)，利用特殊性降低复杂度。
考虑￥中的一个单位圆，我们将其 \(n\) 等分：
用复数 \(\omega_n^k(0\le k<n)\) 表示这些等分线；\(\omega_{n}^k=\omega_n^{k\bmod n}(k\ge n)\)，也就是说 \(\omega_{n}^{k+n}=\omega_n^k(1)\)。这样有什么好处呢？任意两个 \(\omega\) 相乘还在单位圆上，并且：
解释：
\((2)\)：极角相加。
\((3)\)：表示的是同一条线。
\((4)\)：由 \((1)\) 推导而来。事实上，\(\omega_n^k\) 的 \(k\) 表示的就是单位复根（\(w_n^1\)）的 \(k\) 次方。
\((5)\)：类比 \(\alpha\) 和单位圆的交点和 \(\alpha+\pi\) 和单位圆的交点关于原点对称。
这 5 个性质在后面会用到。
FFT——将多项式点值化
将它按下标奇偶性分成两部分：
因此
回顾上述过程，其实是将问题剖分为了两个等规模的子问题——分治。递归地处理 \(A_1,A_2\) 的 FFT 并利用 \((5)\) 归并。
利用「单位复根」相关知识：
当 \(k<\frac n2\) 时
而 \(k+\frac n2\ge \frac n2\)
【注】两式的意义：观察两式，发现只有符号区别，因此是具有统一性的。我们利用了单位复根的特殊性使得问题被成功地分治。
【注2】由于下图的原因，这两个式子统称为“蝴蝶操作”。
IFFT——如何实现点值转系数
引理 将 \(A(x)\) 傅里叶变换的结果 \(A(\omega_{n}^0),A(\omega_n^1),...,A(\omega_n^{n-1})\) 作为一个多项式 \(A’(x)\) 的系数 \(a_0,a_1,...,a_{n-1}\)，代入 \(\omega_n^0,\omega_n^{-1},...,\omega_n^{-(n-1)}\) 傅里叶变换的结果，再除以 \(n\)，得到的就是 \(A(x)\) 的各项系数。
cp omega(int n,int k){
return cp(cos(2*PI*k/n),sin(2*PI*k/n));//分别给z=a+bi的a、b赋值
}
void fft(cp *a,int n,bool inv){
if(n==1)return;
static cp buf[N];
int m=n/2;
for(int i=0;i<m;i++){ //将每一项按照奇偶分为两组
buf[i]=a[2*i];
buf[i+m]=a[2*i+1];
}
for(int i=0;i<n;i++)
a[i]=buf[i];
fft(a,m,inv);//递归处理两个子问题
fft(a+m,m,inv);
for(int i=0;i<m;i++){ //枚举x，计算A(x)
cp x=omega(n,i);
if(inv) x=conj(x);
//conj是一个自带的求共轭复数的函数，精度较高。当复数模为1时，共轭复数等于倒数
buf[i]=a[i]+x*a[i+m];//根据之前推出的结论计算
buf[i+m]=a[i]-x*a[i+m];
}
for(int i=0;i<n;i++)
a[i]=buf[i];
}
//extracted from RabbitHu's Blog
FFT 常数优化
利用下面两个常数优化，有效减少初版 FFT 的运行时间。
非递归 FFT
上面的 FFT 已经是正确的了，但常数较大。递归是罪魁祸首，因此可以从分治后 a 数组的末状态开始倒着算过来，每次还原那一轮的 a。
引理 a[i] 最后会到 a[i’] 去，其中 i’ 为 i 的 \(n\) 位二进制表示（不足补前导0）下翻转的结果，例如 n=3,i=3=011,i’=100=4。
减少数组访问次数
首先，buf 数组的存在就是为了避免调用 a[i],a[i+m] 时会由于前后被修改而互相影响。但如果
cp t=omg[n/l*j]*a[i+j+m];
a[i+j+m]=a[i+j]-t;
a[i+j]+=t;
就可以避免，从而（1）消除了了 buf 数组的无谓使用（2）减少了数组的访问次数。
//https://www.luogu.com.cn/problem/P1919
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N=4e6+5;//第28行的操作会让n最大达到4e6
const double PI=acos(-1); //arccos(-1)=π
typedef complex<double> cp;
char as[N],bs[N];
int n,lgn,la,lb,res[N];
cp omg[N],inv[N],a[N],b[N];
void fft(cp *a,cp *omg){
for(int i=0;i<n;i++){
int t=0;
for(int j=0;j<lgn;j++)t|=((i>>j)&1)<<(lgn-1-j);
if(i<t)swap(a[i],a[t]);
}
for(int l=2;l<=n;l<<=1){
int m=l>>1;
for(int i=0;i<n;i+=l)
for(int j=0;j<m;j++){//omg(n,k)=2*PI*j/n*(n/l)
cp t=omg[n/l*j]*a[i+j+m];
a[i+j+m]=a[i+j]-t;
a[i+j]+=t;
}
}
}
signed main(){
scanf("%s%s",as,bs),la=strlen(as),lb=strlen(bs);
for(n=1;n<la+lb;n<<=1,lgn++);
for(int i=0;i<la;i++)a[i].real(as[la-1-i]-'0'); //.real() 既可以“.real(val)”来给实部赋值，也可以“.real()”返回实部
for(int i=0;i<lb;i++)b[i].real(bs[lb-1-i]-'0');
for(int i=0;i<n;i++)omg[i]=cp(cos(2*PI*i/n),sin(2*PI*i/n)),inv[i]=conj(omg[i]);
fft(a,omg),fft(b,omg);
for(int i=0;i<n;i++)a[i]*=b[i];
fft(a,inv);
for(int i=0;i<n;i++){
res[i]+=a[i].real()/n+0.5; //我一开始不明白这样为什么不会减少精度。
//事实上，理论上这里的实部应该是整数，虚部的系数为0，只是因为复数运算的误差需要四舍五入一下
res[i+1]+=res[i]/10;
res[i]%=10;
}
for(int i=n-1,u=0;i>=0;i--){if(res[i])u=1;if(u)cout<<res[i];}
}
浙公网安备 33010602011771号