震惊！FFT竟然还能这么写？活到这么大没见过这么写FFT的！

模板题

洛谷P3803 【模板】多项式乘法（FFT）

点值表达

大莉模拟的复杂度是$O(n^2)$。
引入点值表达qwq（这部分摘抄自毒瘤czh的FFT讲稿）：
$1.$例子：$A(x)=x^2+2x-1$可以被表达为 { ( 0 , − 1 ) , ( 1 , 2 ) , ( 2 , 7 ) } \{(0 , -1), (1 , 2),(2 , 7)\} {(0,−1),(1,2),(2,7)}

$2.$用处：在多项式乘法中，系数表达暴力模拟复杂度$O(n^2)$，但是点值表达式可以$O(n)$求出

更加具体地说，我们有两个点值表达
A ( x ) = { ( x 0 , A ( x 0 ) ) , ( x 1 , A ( x 1 ) ) , ( x 2 , A ( x 2 ) ) , . . . . , ( x n , A ( x n ) ) } A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),(x_2,A(x_2)),....,(x_n,A(x_n))\} A(x)={(x0​,A(x0​)),(x1​,A(x1​)),(x2​,A(x2​)),....,(xn​,A(xn​))}
B ( x ) = { ( x 0 , B ( x 0 ) ) , ( x 1 , B ( x 1 ) ) , ( x 2 , B ( x 2 ) ) , . . . . , ( x n , B ( x n ) ) } B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),(x_2,B(x_2)),....,(x_n,B(x_n))\} B(x)={(x0​,B(x0​)),(x1​,B(x1​)),(x2​,B(x2​)),....,(xn​,B(xn​))}
乘积 C ( x ) = A ( x ) × B ( x ) = { ( x 0 , A ( x 0 ) × B ( x 0 ) ) , ( x 1 , A ( x 1 ) × B ( x 1 ) ) , ( x 2 , A ( x 2 ) × B ( x 2 ) ) , . . . . , ( x n , A ( x n ) × B ( x n ) ) } C(x)=A(x)\times B(x)=\{(x_0,A(x_0)\times B(x_0)),(x_1,A(x_1)\times B(x_1)),(x_2,A(x_2)\times B(x_2)),....,(x_n,A(x_n)\times B(x_n))\} C(x)=A(x)×B(x)={(x0​,A(x0​)×B(x0​)),(x1​,A(x1​)×B(x1​)),(x2​,A(x2​)×B(x2​)),....,(xn​,A(xn​)×B(xn​))}
所以我们可以通过$O(n)$复杂度，求出两个点值表达式乘积的结果（两个多项式的乘积又称卷积）
知道了点值表达式就珂以确定一个多项式了qwq（珂以列方程用高斯消元爆解一波qwq）

其中DFT算法是暴力$O(n^2)$把系数表达转点值表达，FFT算法是$O(nlogn)$把系数表达转点值表达

DFT（离散傅里叶变换）

普通做法：随机取$n$个值代入求值qwq。
DFT是考虑建立一个复平面，在上面怼出 画一个单位圆：

因为是单位圆，所以圆上的每一个点对应的复数模均为1。
把这个圆平均分成$n$份，把每个点对应的复数代入。
具体地说，这个函数的点值表达是 { ( w n 0 , f ( w n 0 ) ) , ( w n 1 , f ( w n 1 ) ) , . . . , ( w n n − 1 , f ( w n n − 1 ) ) } \{(w_n^0,f(w_n^0)),(w_n^1,f(w_n^1)),...,(w_n^{n-1},f(w_n^{n-1}))\} {(wn0​,f(wn0​)),(wn1​,f(wn1​)),...,(wnn−1​,f(wnn−1​))}，
其中$w_n^k=cos(\frac{k}{n}2\pi )+isin(\frac{k}{n}2\pi )$

FFT（快速傅里叶变换）

$w$的性质

证明珂以看巨神mhy的博客qwq
这里只列出性质，就不证明了：
$1.w_n^k=w_{2n}^{2k}$
$2.w_n^{k+\frac{n}{2}}=-w_n^k$
$3.w_n^0=w_n^n=1$
$4.w_n^p\ast w_n^q=w_n^{p+q}$

表达式的变形

首先把一个函数用系数表达式写出来（假设$n$为偶数，如果不是就补一个$0\ast x^n$）：
$f(x)=a_0+a_1{x}^1+a_2{x}^2+a_3{x}^3+a_4{x}^4+...+a_{n-1}{x}^{n-1}$
按照下标的奇偶性把$f$分成两半：
$f(x)=(a_0+a_2{x}^2+...+a_{n-2}{x}^{n-2})+x(a_1+a_3{x}^2+...+a_{n-1}{x}^{n-2})$
令$A_1(x)=a_0+a_{2}x+a_4{x}^2+...+a_{n-2}{x}^{\frac{n-2}{2}}$，$A_2(x)=a_1+a_{3}x+a_5{x}^2+...+a_{n-1}{x}^{\frac{n-2}{2}}$
不难发现$f(x)=A_1(x^2)+x{A}_2(x^2)$。

以下分别把$w_n^k$和$w_n^{k+\frac{n}{2}}$代入：
把$w_n^k$代入，得到$f(w_n^k)=A_1(w_n^k\ast w_n^k)+x{A}_2(w_n^k\ast w_n^k)=A_1(w_n^{2k})+w_n^k{A}_2(w_n^{2k})$
由$w$的性质1，珂以得到：
$f(w_n^k)=A_1(w_{\frac{n}{2}}^k)+w_n^k{A}_2(w_{\frac{n}{2}}^k)$

把$w_n^{k+\frac{n}{2}}$代入，得到：
$f(w_n^{k+\frac{n}{2}})=A_1(w_n^{2k+n})+w_n^{k+\frac{n}{2}}{A}_2(w_n^{2k+n})=A_1(w_n^{2k}\ast w_n^n)-w_n^k{A}_2(w_n^{2k}\ast w_n^n)=A_1(w_n^{2k})-w_n^k{A}_2(w_n^{2k})$

分治写法

认真观察……这$f(w_n^k)$和$f(w_n^{k+\frac{n}{2}})$长得很像？
于是珂以递归（分治）解决，每次求解左半部分的答案，再推出右半部分的答案：

void FFT(complex<double> *ans,int n,int type) {
	//type=1表示FFT；type=-1表示IFFT 
	if(n==1)	return;
	//偶数存L，奇数存R 
	complex<double> L[(n>>1)+1],R[(n>>1)+1];
	for(re i=0; i<n; i+=2) {
		L[i>>1]=ans[i];
		R[i>>1]=ans[i+1];
	}
	FFT(L,n>>1,type);
	FFT(R,n>>1,type);
	//玄学单位圆乱搞 
	complex<double> w(1,0);
	complex<double> k(cos(2.0*pi/n),type*sin(2.0*pi/n));
	int maxn=n>>1;
	for(re i=0; i<maxn; i++) {
		complex<double> tmp=w*R[i];
		ans[i]=L[i]+tmp;
		ans[i+maxn]=L[i]-tmp;
		w*=k;
	}
}

然后交到例题上，TLE？？
递归常数蜃大qwq，考虑使用迭代。

迭代写法

观察下面这张图：

发现原序列的每个数按对称轴翻转一下就是后序列的数qwq，称之为位逆序置换。
所以，我们可以预处理出系数位置递归到最底层的位置，自底向上进行递推。
处理出位逆序置换之后的方法（$R[i]$表示$i$进行位逆序置换之后得到的值，$L$表示二进制下的总位数）：

for(re i=0; i<MAXN; i++) {		//MAXN表示要处理出来的数的总数
	//表示先把最后一位去掉，剩下的部分进行反转，然后把最后一位反转到最高位qwq
	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
}

所以把递归强行用循环写出来：
第一维枚举区间长度，即递归版本中的$n$。第二维枚举区间开始位置。第三维枚举当前数是区间中的第几个。
然后珂以大莉写一波qwq（我的版本中，区间长度实际上为$2i$）：

void FFT(Complex *ans,int n,int type) {
	//type=1表示FFT，type=-1表示IFFT 
	//把序列交换成后序列的形式 
	for(re i=0; i<n; i++)	if(i<R[i])	swap(ans[i],ans[R[i]]);
	for(re i=1; i<n; i<<=1) {
		//(2*pi)/(2*i)=pi/i
		Complex wn=(Complex){cos(pi/(double)i),type*sin(pi/(double)i)};
		for(re j=0; j<n; j+=(i<<1)) {
			Complex w=(Complex){1,0};
			for(re k=0; k<i; k++) {
				//类似递归版的处理
				const Complex x=ans[j+k],y=w*ans[i+j+k];
				ans[j+k]=x+y;
				ans[i+j+k]=x-y;
				w=w*wn;
			}
		}
	}
}

IFFT（快速傅里叶逆变换）

毒瘤czh：

美妙结论：一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以n即为原多项式的每一项系数（具体看程序实现，证明略）

也就是在调用的时候把type设成-1就珂以了

例题完整代码

至此，就珂以写出例题了（注意预处理和边界条件qwq）：

#include<stdio.h>
#include<cstring>
#include<algorithm>
#include<math.h>
#define re register int
#define rl register ll
#define pi 3.14159265358979323846264338328
using namespace std;
int read() {
    re x=0,f=1;
    char ch=getchar();
    while(ch<'0' || ch>'9') {
        if(ch=='-')    f=-1;
        ch=getchar();
    }
    while(ch>='0' && ch<='9') {
        x=10*x+ch-'0';
        ch=getchar();
    }
    return x*f;
}
inline void write(const int x) {
	if(x>9)	write(x/10);
	putchar(x%10+'0');
}
const int Size=4000005;
//手写complex类，不过貌似也没快多少 
struct Complex {
    double x,y;
    inline double real() {
        return x;
    }
};
inline Complex operator + (const Complex a,const Complex b) {
    return (Complex){a.x+b.x,a.y+b.y};
}
inline Complex operator - (const Complex a,const Complex b) {
    return (Complex){a.x-b.x,a.y-b.y};
}
inline 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};
}
Complex A[Size],B[Size];
int R[Size];
void FFT(Complex *ans,const int n,const int type) {
	//type=1表示FFT，type=-1表示IFFT
	for(re i=0; i<n; i++)	if(i<R[i])	swap(ans[i],ans[R[i]]);
    for(re i=1; i<n; i<<=1) {
    	//(2*pi)/(2*i)=pi/i
        const Complex wn=(Complex){cos(pi/i),type*sin(pi/i)};
        for(re j=0; j<n; j+=i<<1) {
            Complex w=(Complex){1,0};
            for(re k=0; k<i; k++) {
                Complex x=ans[j+k],y=w*ans[i+j+k];
                ans[j+k]=x+y;
                ans[i+j+k]=x-y;
                w=w*wn;
            }
        }
    }
}
int ans[Size];
int main() {
    int n=read();
    int m=read();
    for(re i=0; i<=n; i++) {
        A[i]=(Complex){read(),0};
    }
    for(re i=0; i<=m; i++) {
    	B[i]=(Complex){read(),0};
	}
	//now为大于n+m的第一个2的幂，因为总区间长度必须是2的幂qwq
    int now=1,L=0;
    while(now<=n+m) {
		now<<=1;
		L++;
	}
    for(re i=0; i<now; i++) {
    	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
	}
    FFT(A,now,1);
    FFT(B,now,1);
    for(re i=0; i<=now; i++) {
        A[i]=A[i]*B[i];
    }
    FFT(A,now,-1);
    for(re i=0; i<=n+m; i++) {
		write(int(A[i].real()/now+0.5));	//注意四舍五入
		putchar(' ');
    }
    return 0;
}

其他题目

洛谷P1919 A*B Problem升级版
即bzoj2179 （权限题）快速傅里叶变换

bzoj2194（权限题）快速傅里叶变换之二

bzoj3527 力

bzoj4827 礼物
……
…………
………………
题解？咕咕咕

还是有一篇的qwq

bzoj4332（权限题）
即洛谷P5075
题解传送门