[学习笔记]快速傅里叶变换
一. 多项式
1. 定义
我们称一个关于 \(x\) 的式子
为一个 \(n\) 次多项式，其中 \(a_i\) 称为系数。
称 \(n\) 为 \(f(x)\) 的次数，那么显然地 \(f(x)\) 可以看做一个 \(x\) 的一元 \(n\) 次方程
那么比较显然地，想确定一个 \(n\) 次函数的解析式需要至少 \(n+1\) 个点的坐标。
可以用这 \(n+1\) 个点表示 \(f(x)\) ，称为点值表达法
对应地，把上面的那个式子称作系数表示法
2. 基本运算规律
以下我们设两个多项式 \(f(x)=\sum\limits_{i=0}^{n}a_ix^i\) ，\(g(x)=\sum\limits_{i=0}^{n}a_ix^i\)
\(\qquad\) \(\mathrm{I.}\) 加减法：
\(\qquad\qquad\) 显然地，\(f(x)\pm g(x)=\sum\limits_{i=0}^{n}(a_i+b_i)x^i\)
\(\qquad\) \(\mathrm{II.}\) 乘法：
\(\qquad\qquad\) 两个 \(n\) 次多项式的乘积是一个 \(2n\) 次的多项式
\(\qquad\qquad\) 也可以写作
二. 复数
请见高中数学必修二第七章《复数》，这里只提供简要讲解。
1. 定义
定义常数 \(i\) ，满足 \(i^2=-1\)
则所有形如 \(z = a + b \times i\;\; (a, b \in \mathfrak{R})\) 的数字 \(z\) 构成的集合称为复数集，记为 \(\mathbb{C}\)
\(\mathbb{C}\) 中的每个元素都称作复数。
对于 \(z=a+bi\) ，称 \(a\) 是 \(z\) 的实部，\(b\) 是 \(z\) 的虚部。
2. 几何意义
复平面是一个笛卡尔平面，有两条坐标轴，纵轴为虚轴，横轴为实轴，两轴相互垂直。
对于一个复数 \(z=a+bi\) ，它在复平面上对应一个从原点指向 \((a,b)\) 的向量。也即实轴坐标为实部值，虚轴坐标为虚部值的点。
显然复平面上从原点出发的任意一个向量也对应唯一的一个复数，也即向量 \((a,b)\) 对应一个复数 \(z=a+bi\)
易证复数与复平面上从原点出发的向量是一一对应关系。
以实轴正方向为始边，\(z\) 所对应的向量 \(\vec{z}\) 为终边的角 \(\theta\) 称为复数 \(z\) 的幅角。
对于所有的幅角，我们把在 \([0,2\pi]\) 间的幅角称为 \(z\) 的幅角主值，记作 \(\arg z\)
这样我们还可以用三角函数表示一个复数 \(z\) ：\(z=\left\vert z\right\vert(\cos(\arg z)+i\sin(\arg(z)))\)
3.运算律
\(\qquad\) \(\mathrm{I.}\) 复数的模：
\(\qquad\qquad\) 复数 \(z=a+bi\) 的模为其在复平面上对应向量的长度，记作 \(\left\vert z\right\vert\)
\(\qquad\qquad\) \(\Large\left\vert z\right\vert=\sqrt{a^2+b^2}\)
\(\qquad\) \(\mathrm{II.}\) 共轭复数：
\(\qquad\qquad\) 复数 \(z\) 在复平面上对应的向量关于实轴对称后对应的复数称为 \(z\) 的共轭复数，记作 \(\overline{z}\)
\(\qquad\qquad\) \(z\) 与 \(\overline{z}\) 满足如下关系：
\(\qquad\qquad\qquad\) \(①\) 设 \(\arg z=\theta_0\) ，\(\arg\overline{z}=\theta_1\) ，则 \(\theta_0 + \theta_1 = \pi\)
\(\qquad\qquad\qquad\) \(②\) \(\left\vert z\right\vert = \left\vert\overline{z}\right\vert\)
\(\qquad\qquad\qquad\) \(③\) \(z\) 和 \(\overline{z}\) 的实部相同，虚部互为相反数。
\(\qquad\) \(\mathrm{III.}\) 加减法：
\(\qquad\qquad\) 不妨设 \(z_1=a_1+b_1i\) ，\(z_2=a_2+b_2i\)
\(\qquad\qquad\) 显然地，\(z_1 \pm z_2~=~(a_1 \pm a_2) + (b_1 \pm b_2)i\)
\(\qquad\qquad\) 在复平面上，两个复数相加减的结果为他们所对应的向量按照平行四边形定则相加减。
\(\qquad\) \(\mathrm{IV.}\) 乘法：
\(\qquad\qquad\) \(z_1=a_1+b_1i\) ，\(z_2=a_2+b_2i\)
\(\qquad\qquad\) \(\Large\begin{aligned} z_1 \times z_2 & =(a_1 + b_1i) \times (a_2 + b_2i) = a_1a_2 + a_1b_2i + a_2b_1i + b_1b_2i^2\\ & =(a_1a_2 - b_1b_2) + (a_1b_2 + a_2b_1) i\end{aligned}\)
\(\qquad\qquad\) 在复平面上，\(z_1,z_2,z_0=z_1\times z_2\) 所对应的幅角 \(\theta_1,\theta_2,\theta_0\) 有如下关系：
\(\qquad\qquad\qquad\) \(\large\theta_0 = \theta_1 + \theta_2\)
\(\qquad\qquad\) 他们的模有如下关系：
\(\qquad\qquad\qquad\) \(\large\left\vert z_0\right\vert= \left\vert z_1\right\vert \times \left\vert z_2\right\vert\)
\(\qquad\qquad\) 也就是幅角相加，模长相乘
\(\qquad\qquad\) 考虑 \(z = a + b_i\) 和它的共轭复数 \(\overline{z} = a - bi\)
\(\qquad\qquad\) \(\large z \times \overline{z}=(a + bi) \times (a - bi) = a^2 + b^2\)
\(\qquad\qquad\) 因此两个互为共轭复数的数之积一定是一个实数。
\(\qquad\) \(\mathrm{V.}\) 除法：
\(\qquad\qquad\) 对于两个复数 \(z_1 = a_1 + b_1i,z_2 = a_2 + b_2 i\)，他们相除的结果为 \(\frac{z_1}{z_2}\)
\(\qquad\qquad\) 考虑分数上下同时乘 \(\overline{z_2}\) ，有 \(\Large\frac{z_1}{z_2}=\frac{z_1 \overline {z_2}}{a_2^2 + b_2^2}\)
\(\qquad\qquad\) 分母是一个实数，可以直接将分子的实部虚部除以分母。
\(\qquad\) \(\mathrm{VI.}\) 复数指数幂：
\(\qquad\qquad\) 有欧拉公式
\(\qquad\qquad\) 其中 \(e\) 是自然对数的底数（也即 \(e=\lim\limits_{x\to +\infty}(1+x)^\frac{1}{x}）\)
三. 单位根
在复数域内，满足 \(x^n=1\) 的 \(x\) 称为 \(x\) 的 \(n\) 次单位根。
根据代数基本定理， \(n\) 次单位根一共有 \(n\) 个。
经过计算可得，将所有的 \(n\) 次单位根按照幅角大小排列，第 \(k\;(0 \leqslant k < n)\) 个 \(n\) 次单位根为
因为 \(\sin^2 \theta + \cos^2 \theta=1\) ，所以所有的 \(n\) 次单位根的模都是 \(1\)
\(n\) 个单位根在复平面上平分单位圆，他们与单位圆的 \(n\) 等分线所在线段分别重合。
四. 本原单位根
对于一个 \(n\) 次单位根，如果它的 \(0\) 到 \((n−1)\) 次方的值能生成所有 \(n\) 次单位根，就把它称为 \(n\) 次本原单位根。
显然地，\(\large x_1 = e^{i \frac{2\pi}{n}}\) 是一个本原单位根
我们把 \(n\) 次本原单位根称为 \(\large \omega_n=e^{i \frac{2\pi}{n}}=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}\)
那么不妨设 \(n=2m\)，则有
五. 离散傅里叶变换 (DFT)
首先根据小学数学，我们知道 \(n+1\) 个点能唯一确定一个 \(n\) 次多项式
假设 \(f(x)\) 和 \(g(x)\) 都是 \(n\) 次多项式
给了你 \(f(x)\) 的 \(n+1\) 个点 \((x_1,y_{f1})~,~(x_2,y_{f2})~,~\cdots~,~(x_n,y_{fn})\) 和 \(g(x)\) 的 \(n+1\) 个点 \((x_1,y_{g1})~,~(x_2,y_{g2})~,~\cdots~,~(x_n,y_{gn})\)
考虑怎样求出 \(h(x)=f(x)\times g(x)\)
根据定义，\(h(x)\) 在 \(x=x_i\) 处的值就是 \(y_{fi}\times y_{gi}\)
很好，现在我们就求出了 \(h(x)\)
但是显然的是，一般人并不会用 \(n+1\) 个点来表示一个多项式
我们想要一个用系数表示的式子
所以我们考虑
六. 快速傅里叶变换 (FFT)
首先我们设 \(n\) 为 \(2\) 的整数次幂
等等，你说有 \(n\) 不为 \(2\) 的整数次幂的情况？
那我们就补 \(0\) 到合法为止
再设一个 \(\large m=\frac{n}{2}\)
设 \(f(x)\) 为一个 \(n-1\) 次多项式，注意到
不妨设 \(\large f_1(x)=\sum\limits_{i=0}^{m-1}a_{2i}x^{i}~,~f_2(x)=\sum\limits_{i=0}^{m-1}a_{2i+1}x^{i}\)
那么显然地，\(\large f(x)=f_1(x^2)+xf_2(x^2)\)
不妨代入 \(\omega_n^k\; (k<m)\)
根据上面单位根的性质，我们有
再代入 \(\omega_n^{k+m}\; (k<m)\)
也就是说，理论上来讲，只要我们知道 \(f_1(x)\) 和 \(f_2(x)\) 在 \(\omega_m^0~,~\omega_m^1~,~\cdots~,~\omega_m^{m-1}\) 处的点值表示
就可以 \(\mathcal{O(n)}\) 求出 \(f(x)\) 的点值表达式
然后我们惊讶地发现：原来 \(f_1(x)\) 和 \(f_2(x)\) 是可以像 \(f(x)\) 一样处理的？
递归就好了
根据主定理，时间复杂度是 \(\mathcal{O}(n\log n)\)
上面这玩意就是 \(\operatorname{FFT}\) 狒狒贴
于是我们可以写出一个叫做 \(\operatorname{complex}\) 的结构体用来存储复数
复数的四则运算
struct Complex{
double x,y;
Complex(double x_,double y_){x=x_,y=y_;}
Complex(){}
double magnitude(){//求模长
return x*x+y*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);
}
Complex operator/(const Complex &b)const{
double mag=b.magnitude();
return Complex((x*b.x+y*b.y)/mag,y*b.x-x*b.y);
}
};
然后还可以顺便求一下单位根们：
求单位根
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define int long long
#define WR WinterRain
using namespace std;
const int WR=2001000;
const double pi=acos(-1);
struct Complex{
double x,y;
Complex(double x_,double y_){x=x_,y=y_;}
Complex(){}
double magnitude(){//求模长
return x*x+y*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);
}
Complex operator/(Complex b){
double mag=b.magnitude();
return Complex((x*b.x+y*b.y)/mag,y*b.x-x*b.y);
}
}omega[WR];
int n;
int read(){
int s=0,w=1;
char ch=getchar();
while(ch>'9'||ch<'0'){
if(ch=='-') w=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
s=(s<<3)+(s<<1)+ch-'0';
ch=getchar();
}
return s*w;
}
signed main(){
n=read();
Complex base=Complex(cos(2.0*pi/n),sin(2.0*pi/n));
omega[0]=Complex(1,0);
for(int i=1;i<n;i++){
omega[i]=omega[i-1]*base;
}
for(int i=0;i<n;i++){
printf("omega[%lld]=(%.5lf,%.5lfi)\n",i,omega[i].x,omega[i].y);
}
return 0;
}
然后就可以进行一个 \(\operatorname{DFT}\) 的写
DFT（递归版）
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define int long long
#define WR WinterRain
using namespace std;
const int WR=2001000;
const double pi=acos(-1);
struct Complex{
double x,y;
Complex(double x_,double y_){x=x_,y=y_;}
Complex(){}
double magnitude(){//求模长
return x*x+y*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);
}
Complex operator/(Complex b){
double mag=b.magnitude();
return Complex((x*b.x+y*b.y)/mag,y*b.x-x*b.y);
}
}rec[WR];
int n,lenth;
int read(){
int s=0,w=1;
char ch=getchar();
while(ch>'9'||ch<'0'){
if(ch=='-') w=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
s=(s<<3)+(s<<1)+ch-'0';
ch=getchar();
}
return s*w;
}
void DFT(Complex *f,int len){
if(len==1) return;
int m=(len>>1);
Complex *fl=f,*fr=f+m;
for(int i=0;i<len;i++) rec[i]=f[i];
for(int i=0;i<m;i++){
fl[i]=rec[i<<1];
fr[i]=rec[i<<1|1];
}
DFT(fl,len>>1);//分治
DFT(fr,len>>1);
Complex omega=Complex(cos(2.0*pi/len),sin(2.0*pi/len)),base=Complex(1,0);
for(int i=0;i<m;i++){
rec[i]=fl[i]+base*fr[i];
rec[i+m]=fl[i]-base*fr[i];
base=base*omega;
}
for(int i=0;i<len;i++) f[i]=rec[i];
}
Complex f[WR];
signed main(){
n=read();n++;
for(int i=0;i<n;i++) scanf("%lf",&f[i].x);
for(lenth=1;lenth<n;lenth<<=1);//找到 2 的整数次幂
DFT(f,lenth);
for(int i=0;i<lenth;i++) printf("(%.5lf,%.5lf)\n",f[i].x,f[i].y);
return 0;
}
上面的 \(\operatorname{DFT}\) 成功地将一个系数表示的多项式转化成了点值表达
但遗憾的是，我们要求的是系数表达
怎么办呢？
考虑下面的内容
七. 傅里叶逆变换（ \(\operatorname{IDFT}\) ）
首先有一个叫做范德蒙矩阵的玩意长这样
它有一些美妙的性质
设 \(\large f(x)=\sum\limits_{i=0}^{n}a_ix^i\) ，我们有一个系数矩阵
根据矩阵乘法原理，显然地，
也就是说，用一个系数矩阵乘上范德蒙矩阵相当于多点求值
相应地，用一个多项式值矩阵乘上范德蒙矩阵的逆矩阵相当于多点插值
不妨设多项式 \(\large\sum\limits_{i=0}^{n-1}a_ix^i\)
不难发现，用矩阵表示的 \(\operatorname{DFT}\) 是：
相应地，我们要求的是
现在我们考虑证明这个范德蒙矩阵 \(V(\omega_n^0~,~\omega_n^1~,~\omega_n^{n-1})\) 是可求逆的（也即 \(\left\vert V\right\vert\neq 0\) )
结论：范德蒙行列式 \(V(x_0~,~x_1~,~\cdots~,~x_n)\) 的值为 \(\left\vert V\right\vert=\prod\limits_{1\leqslant j<i\leqslant n}(x_i-x_j)\; (~n\geqslant 2~)\)
采用数学归纳法
首先对于 \(n=2\) ，\(V(x_0~,~x_1)=\begin{vmatrix} 1 & 1\\x_0 & x_1\\ \end{vmatrix}=x_1-x_0\) 显然成立
在该结论对于 \(n\) 阶以下（包括 \(n\) 阶）的范德蒙矩阵全部成立的时候，
首先注意到有一个 \(n\) 阶范德蒙矩阵
不妨设其为 \(V_{pre}\)
然后我们考察 \(V(x_0~,~x_1~,~\cdots~,~x_n)\)
考虑对其运用初等行变换
按照第一列展开，发现消去了绝大多数内容，只留下
继续初等行变换，提取公因式得出
最后，由于显然地任意两个 \(n\) 次单位根不可能相等，也就是说这个行列式是非奇异的
可以求逆。如何求逆？
不妨设上面的范德蒙行列式为 \(V\)
考虑如何多点插值，自然地想到拉格朗日插值
注意到有拉格朗日插值公式
容易发现，\(V^{-1}_{i,j}\) 就是 \(\prod_{k\neq i}\frac{x-x_k}{x_i-x_k}\) 在 \(x_{j-1}\) 这项的系数
现在我们有
不妨令 \(g(x)=\prod\limits_{0\leqslant j<n}(x-\omega_n^j)\)
而注意到 \(\omega_n^i\) 都是 \(1\) 的 \(n\) 次根
根据代数基本定理，\(g(x)=x^n-1\)
分离后面的连乘，可以化成
分别考察分子分母：
对于分母，显然是一个 \(\frac{0}{0}\) 的未定式
考虑洛必达法则：
对于分子，注意到
考虑对其泰勒展开，我们有泰勒公式
现在令 \(f(x)=\frac{1}{1-x}\)
显然地，
\(f^\prime(x)=\left(\frac{1}{x}\right)^\prime(1-x)^\prime=\frac{1!}{(1-x)^2}\)
\(f^{\prime\prime}(x)=\left(\frac{1}{(1-x)^2}\right)^\prime=\frac{2!}{(1-x)^3}\)
\(\cdots\cdots\)
\(f^{(n)}(x)=\frac{n!}{(1-x)^{n+1}}\)
将 \(x_0=0\) 代入，可得
代回上式，得到
考察变化后的分子分母，我们可以得到
好眼熟啊
等等，这不就是上面 \(\operatorname{DFT}\) 的单位根求个了倒数吗？
然后再乘一个 \(\frac{1}{n}\) 即可
也就是说，\(\operatorname{IDFT}\) 和 \(\operatorname{DFT}\) 只有一个正负号的区别！
八. FFT代码实现
显然可以写出
点击查看代码
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define int long long
#define WR WinterRain
using namespace std;
const int WR=5001000;
const double pi=acos(-1);
struct Complex{
double x,y;
Complex(double x_,double y_){x=x_,y=y_;}
Complex(){}
double magnitude(){
return x*x+y*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);
}
Complex operator/(Complex b){
double mag=b.magnitude();
return Complex((x*b.x+y*b.y)/mag,y*b.x-x*b.y);
}
}rec[WR];
int read(){
int s=0,w=1;
char ch=getchar();
while(ch>'9'||ch<'0'){
if(ch=='-') w=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
s=(s<<3)+(s<<1)+ch-'0';
ch=getchar();
}
return s*w;
}
void DFT(Complex *f,int len){
if(len==1) return;
int m=(len>>1);
Complex *fl=f,*fr=f+m;
for(int i=0;i<len;i++) rec[i]=f[i];
for(int i=0;i<m;i++){
fl[i]=rec[i<<1];
fr[i]=rec[i<<1|1];
}
DFT(fl,len>>1);
DFT(fr,len>>1);
Complex omega=Complex(cos(2.0*pi/len),sin(2.0*pi/len)),base=Complex(1,0);
for(int i=0;i<m;i++){
rec[i]=fl[i]+base*fr[i];
rec[i+m]=fl[i]-base*fr[i];
base=base*omega;
}
for(int i=0;i<len;i++) f[i]=rec[i];
}
void IDFT(Complex *f,int len){
if(len==1) return;
int m=(len>>1);
Complex *fl=f,*fr=f+m;
for(int i=0;i<len;i++) rec[i]=f[i];
for(int i=0;i<m;i++){
fl[i]=rec[i<<1];
fr[i]=rec[i<<1|1];
}
IDFT(fl,len>>1);
IDFT(fr,len>>1);
Complex omega=Complex(cos(2.0*pi/len),-sin(2.0*pi/len)),base=Complex(1,0);//注意这里有一个负号
for(int i=0;i<m;i++){
rec[i]=fl[i]+base*fr[i];
rec[i+m]=fl[i]-base*fr[i];
base=base*omega;
}
for(int i=0;i<len;i++) f[i]=rec[i];
}
int n,m,len=1;
Complex f[WR],g[WR];
signed main(){
n=read(),m=read();
for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
while(len<=n+m) len<<=1;
DFT(f,len),DFT(g,len);
for(int i=0;i<len;i++) f[i]=f[i]*g[i];
IDFT(f,len);
for(int i=0;i<=n+m;i++) printf("%lld ",(int)(f[i].x/len+0.49));
printf("\n");
return 0;
}
但是可以看到时间 \(2.5s\)
能不能快一点啊？
首先压个行，把 \(\operatorname{DFT}\) 和 \(\operatorname{IDFT}\) 合并
然后注意到了我们用指针进行数组的拷贝
这无疑是非常浪费时间的，考虑进行优化
不难发现最后序列的位置是原序列的 二进制反转
也就是说原来的 \(6=(110)_2\) 在变换到最底层时在 \((011)_2=3\) 的位置上
不妨设一个二进制数为 \((\overline{abc\cdots x})_2\)
再开一个数组 \(trans\) 存反转完的数
翻转的话就相当于 \(x\) 连接上其他部分的反转
其他部分的反转就是 \(r[i>>1]>>1\)
然后判一下最后一位,如果是 \(1\) 则加上 \(n>>1\)
for(int i=0;i<len;i++){
trans[i]=(trans[i>>1]>>1)|((i&1)?(len>>1):0);
}
然后注意到有递归，是坏的
考虑对其使用迭代实现
点击查看代码
#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define int long long
#define WR WinterRain
using namespace std;
const int WR=5001000;
const double pi=acos(-1);
struct Complex{
double x,y;
Complex(double x_,double y_){x=x_,y=y_;}
Complex(){}
double magnitude(){
return x*x+y*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);
}
Complex operator/(Complex b){
double mag=b.magnitude();
return Complex((x*b.x+y*b.y)/mag,y*b.x-x*b.y);
}
}rec[WR];
int trans[WR];
int read(){
int s=0,w=1;
char ch=getchar();
while(ch>'9'||ch<'0'){
if(ch=='-') w=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9'){
s=(s<<3)+(s<<1)+ch-'0';
ch=getchar();
}
return s*w;
}
void FFT(Complex *f,int len,bool tag){
for(int i=0;i<len;i++){
if(i<trans[i]) swap(f[i],f[trans[i]]);
}
for(int i=2;i<=len;i<<=1){//枚举区间长度
int m=i>>1;//枚举待合并的长度
Complex omega=Complex(cos(2.0*pi/i),sin(2.0*pi/i));
if(!tag) omega.y*=-1.0;
for(int j=0;j<len;j+=i){//枚举起始点
Complex base=Complex(1.0,0.0);
for(int k=j;k<j+m;k++){//处理左半部分
Complex tmp=base*f[k+m];
f[k+m]=f[k]-tmp;
f[k]=f[k]+tmp;
base=base*omega;
}
}
}
}
int n,m,len=1;
Complex f[WR],g[WR];
signed main(){
n=read(),m=read();
for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
for(int i=0;i<=m;i++) scanf("%lf",&g[i].x);
while(len<=n+m) len<<=1;
for(int i=0;i<len;i++){
trans[i]=(trans[i>>1]>>1)|((i&1)?(len>>1):0);
}
FFT(f,len,1),FFT(g,len,1);
for(int i=0;i<len;i++) f[i]=f[i]*g[i];
FFT(f,len,0);
for(int i=0;i<=n+m;i++) printf("%lld ",(int)(f[i].x/len+0.49));
printf("\n");
return 0;
}
发现效果拔群，整整快了一倍！
这样就差不多是 \(\operatorname{FFT}\) 的简要内容了，希望对你有帮助qwq
本文来自博客园，作者：冬天丶的雨，转载请注明原文链接：https://www.cnblogs.com/WintersRain/p/16685902.html
为了一切不改变的理想，为了改变不理想的一切