FFT学习笔记

最近终于能干自己想干的事情了

问题引入:

快速求多项式乘法（卷积）

引理:

n个不同点的点集可以确定唯一 的n次多项式

记相乘的两个多项式为$A,B$
$A,B$的卷积记为$C$

solution1

$O(n^2)$的多项式乘法

solution2

根据引理,可以任意取$A,B$上 $|A|+|B|$个点
然后把 其纵坐标相乘，故在通过这$|A|+|B|$个点来求出C
复杂度 $O(n^3)[高斯消去]$ / $O(n^2)[拉格朗日插值]$/…

solution3

瓶颈:通过点值表达式来求出$C$
优化:每次可以选择一些横坐标不同，可以用速算点加速
复杂度$O(nlogn)$
实现具体如下:

重要概念

令$W(n,i)$为将单位圆平分成n分的第i个点 $(0<=i<n)$
$W(2\ast n,\mathrm{２}＊i)＝W(n,i)$
$W(n,i+j)＝W(n,i)\ast W(n,j)$
$W(n,i+n/2)=-W(n,i)$

重要公式

$e^{ix}=cosx+isinx$

联系$W(n,i)$，可以得到:

$W(n,i)=cos(2\ast PI\ast i/n)+isin(2\ast PI\ast i/n)$

正戏开始

$A(w_n^k)={\sum }_{i=0}^{n-1}{a}_i\ast W_n^{k+i})$
$B(w_n^k)={\sum }_{i=0}^{\frac{n-1}{2}}{a}_{i\ast 2}\ast W_n^{k\ast i\ast 2})$
$C(w_n^k)={\sum }_{i=0}^{\frac{n-1}{2}}{a}_{i\ast 2+1}\ast W_n^{k(i\ast 2+1)})=$
$W_n^k\ast {\sum }_{i=0}^{\frac{n-1}{2}}{a}_{i\ast 2+1}\ast W_n^{k\ast i\ast 2})$
$A(w_n^k)=B(w_n^k)+C(W_n^k)$

同理

$A(w_n^{k+\frac{n}{2}})={\sum }_{i=0}^{n-1}{a}_i\ast W_n^{k+i})$
$B(w_n^{k+\frac{n}{2}})={\sum }_{i=0}^{\frac{n-1}{2}}{a}_{i\ast 2}\ast W_n^{(k+\frac{n}{2})\ast i\ast 2})$
$C(w_n^{k+\frac{n}{2}})={\sum }_{i=0}^{\frac{n-1}{2}}{a}_{i\ast 2+1}\ast W_n^{k(k+\frac{n}{2})(i\ast 2+1)})=$
$(-1)\ast W_n^k\ast {\sum }_{i=0}^{\frac{n-1}{2}}{a}_{i\ast 2+1}\ast W_n^{(k+\frac{n}{2})+i\ast 2})$

综上

$A(x)=B(x)+C(x)$
$A(w_n^k)=B(w_n^{2k})+C(w_n^{2k})$
$A(w_n^{k+\frac{n}{2}})=B(w_n^{k+\frac{n}{2}})+C(w_n^{k+\frac{n}{2}})$

人话:

就是根据$w_n^i$的循环性质，一次求出每个$x$取$w_n^i$的时候整个多项式的值，求出点值表达式，就在$O(nlogn)$的时间解决了这个问题

点值表达式转回原表达式

直接每次把$W_n^i$取倒，而后做一次fft，且每一项都$\ast n^{-1}$

具体如下

#include<bits/stdc++.h>
#define MAXN 4000005
using namespace std;

int n,m,len,limit,rev[MAXN];
const double PI = acos(-1.0);
struct node{double shi,xu;}a[MAXN],b[MAXN];
node operator + (node x , node y){return (node){x.shi + y.shi , x.xu + y.xu};}
node operator - (node x , node y){return (node){x.shi - y.shi , x.xu - y.xu};}
node operator * (node x , node y){return (node){x.shi * y.shi - x.xu * y.xu , x.shi * y.xu + x.xu * y.shi};}

node dw(double x , double y , int h){//W(x , y) 
	double zz = y / x;
	zz = 2.0 * zz * PI;
	return (node){cos(zz) , sin(zz) * h};
}

void fft(node f[] , int tp){
	for(int i = 0 ; i < len ; i++)if(i < rev[i])swap(f[i] , f[rev[i]]);
	for(int i = 1 ; i < len ; i = i * 2){
		for(int j = 0 ; j < len ; j = j + (i * 2)){
			node xs = dw(i * 2 , 1 , tp) , W = (node){1 , 0};
			for(int k = 0 ; k < i ; k++ , W = W * xs){
				node p = f[k + j] , pp = f[k + j + i] * W;
				f[k + j] = p + pp;
				f[k + j + i] = p - pp;
			}
		}
	}
}

int main(){
	cin>>n>>m , n++ , m++;
	for(int i = 0 ; i < n ; i++)cin>>a[i].shi;
	for(int i = 0 ; i < m ; i++)cin>>b[i].shi;
	while((1 << limit) <= n + m)limit++;
	len = (1 << limit);
	for(int i = 0 ; i < len ; i++){
		rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (limit - 1)));
	}
	fft(a , 1) , fft(b , 1);
	for(int i = 0 ; i < len ; i++)a[i] = a[i] * b[i];
	fft(a , -1);
	for(int i = 0 ; i <= (n + m - 2) ; i++)cout<<int((a[i].shi / (len * 1.0) + 0.5))<<" ";
}