数学板块学习之FFT

FFT

FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。—baidu

信息学竞赛中FFT用于快速求解多项式乘法，也就是卷积

定义两个向量a,b,其中$a=(a_0,a_1,\cdots ,a_{n-1}),b=(b_0,b_1,\cdots ,b_{n-1})$
向量和$a+b=(a_0+b_0,a_1+b_1,\cdots ,a_{n-1}+b_{n-1})$
点积$a\cdot b=(a_0{b}_0,a_1{b}_1,\cdots ,a_{n-1}{b}_{n-1})$
卷积$a\ast b=(c_0,c_1,\cdots ,c_{2n-2}),c_k={\sum }_{i+j=k}{a}_i\cdot =b_j$

系数表示法与点值表示法
给出一个多项式$A(x)$
系数表示法：$A(x)=a_{0}x+a_1{x}^1+\cdots +a_{n-1}{x}^{n-1}={\sum }_{i=0}^{n-1}{a}_i{x}^i\iff$ 系数向量$a=(a_0,a_1,\cdots ,a_{n-1})$
点值表示法： A ( x ) = { ( x 0 , y 0 ) , ( x 1 , y 1 ) , ⋯   , ( x n − 1 , y n − 1 ) } ( x 0 , x 1 , ⋯   , x n − 1 ) A(x)=\{(x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\}(x_0,x_1,\cdots,x_{n-1}) A(x)={(x0​,y0​),(x1​,y1​),⋯,(xn−1​,yn−1​)}(x0​,x1​,⋯,xn−1​)($x_i$互不相同)

转换
系数表示法–>点值表示法
$$y_i=\sum _{k=0}^{n-1}{c}_i{x}_i^k$$
点值表示法–>系数表示法
$$f(x)=\sum _{i=0}^{n-1}{y}_i\frac{{\Pi }_{j≠i}^{0\le j<n}(x-x_j)}{{\Pi }_{j≠i}^{0\le j<n}(x_i-x_j)}$$
由点值转换为系数是拉格朗日插值法

如果多项式$A(x),B(x)$(假设两个多项式的次数和小于$n$,即$degA+degB<n$)的点值表示都是在 { x 0 , x 1 , ⋯ &ThinSpace; , x n − 1 } \{x_0,x_1,\cdots,x_{n-1}\} {x0​,x1​,⋯,xn−1​}处
A ( x ) = { ( x 0 , A ( x 0 ) ) , ( x 1 , A ( x 1 ) ) , ⋯ &ThinSpace; , ( x n − 1 , A ( x n − 1 ) ) } A(x)=\{(x_0,A(x_0)),(x_1,A(x_1)),\cdots,(x_{n-1},A(x_{n-1}))\} A(x)={(x0​,A(x0​)),(x1​,A(x1​)),⋯,(xn−1​,A(xn−1​))}
B ( x ) = { ( x 0 , B ( x 0 ) ) , ( x 1 , B ( x 1 ) ) , ⋯ &ThinSpace; , ( x n − 1 , B ( x n − 1 ) ) } B(x)=\{(x_0,B(x_0)),(x_1,B(x_1)),\cdots,(x_{n-1},B(x_{n-1}))\} B(x)={(x0​,B(x0​)),(x1​,B(x1​)),⋯,(xn−1​,B(xn−1​))}
( A ⋅ B ) ( x ) = { ( x 0 , A ( x 0 ) ⋅ B ( x 0 ) ) , ( x 1 , A ( x 1 ) ⋅ B ( x 1 ) ) , ⋯ &ThinSpace; , ( x n − 1 , A ( x n − 1 ) ⋅ B ( x n − 1 ) ) } (A·B)(x)=\{(x_0,A(x_0)·B(x_0)),(x_1,A(x_1)·B(x_1)),\cdots,(x_{n-1},A(x_{n-1})·B(x_{n-1}))\} (A⋅B)(x)={(x0​,A(x0​)⋅B(x0​)),(x1​,A(x1​)⋅B(x1​)),⋯,(xn−1​,A(xn−1​)⋅B(xn−1​))}
即$(A\cdot B)(x_i)=A(x_i)\cdot B(x_i)$

我们可以通过n+1个点确定两个n阶多项式是否相等:
定理：如果$f(x),g(x)$是$R[x]$的两个多项式，且两个多项式的次数都不大于n,若以$R$中的n+1个点或更多的点代替$x$得到的每个$f(x),g(x)$都相等，则$f(x)=g(x)$

由该定理知道我们可以通过n个点确定一个n-1阶多项式

引入单位复根
我们知道世界上最完美的公式就是欧拉公式$e^{i\pi }=-1$
对于任意实数有$e^{ix}=cos(x)+isin(x)$

我们称满足${\omega }^n=1$的$\omega$为1的n次复根，其中${\omega }_n^k=e^{\frac{ik2\pi }{n}}$,k为从x轴正半轴逆时针旋转的编号，明显k为n时${\omega }_n^n=1$，k为$\frac{n}{2}$时${\omega }_n^{\frac{n}{2}}=-1$
因此对于${\omega }_n^k$我们只需要求出${\omega }_n^1$就可以得到${\omega }_n^k$，所以我们称${\omega }_n^1$为单位复根，简写为${\omega }_n$
性质：
${\omega }_{2n}^{2k}={\omega }_n^k$
${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k$

我们将单位根${\omega }_n$代替x带入多项式，也就是带入${\omega }_n^0,{\omega }_n^1,\cdots ,{\omega }_n^{n-1}$，得到一种特殊的点值表示，这就是傅里叶变换
即$$y_k=\sum _{i=0}^{n-1}{c}_i({\omega }_n^k{)}^i$$
$$\left(\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ⋮\\ y_{n-1}\end{array}\right)=\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_n^1& {\omega }_n^2& \cdots & {\omega }_n^{n-1}\\ 1& {\omega }_n^2& {\omega }_n^4& \cdots & {\omega }_n^{2(n-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }_n^{n-1}& {\omega }_n^{2(n-1)}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right)\left(\begin{array}{c}{c}_0\\ c_1\\ c_2\\ ⋮\\ c_{n-1}\end{array}\right)$$
而逆变换$$\begin{array}{rl}{c}_k& =\frac{1}{n}\sum _{i=0}^{n-1}{y}_i({\omega }_n^{-k}{)}^i\\ & =\frac{1}{n}\sum _{i=0}^{n-1}(\sum _{j=0}^{n-1}{c}_j({\omega }_n^i{)}^j)({\omega }_n^{-k}{)}^i\\ & =\frac{1}{n}\sum _{j=0}^{n-1}{c}_j(\sum _{i=0}^{n-1}({\omega }_n^{j-k}{)}^i)\\ & =c_k\end{array}$$
故傅里叶逆变换只需要改变傅里叶变换的参数$1$变为$-1$并除以$\frac{1}{n}$即可
FFT

说了这么多前置，终于开始了FFT
FFT采用分治的思想，将奇数项和偶数项分开。
对于一个多项式
$$A(x)=(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+(a_1{x}^1+a^3{x}^3+\cdots +a_{n-1}{x}^{n-1})$$
则
$$A(x)=(a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2})+x(a_1+a^3{x}^2+\cdots +a_{n-1}{x}^{n-2})$$
我们设$A_1(x)=a_0+a_2{x}^2+\cdots +a_{n-2}{x}^{n-2},A_2(x)=(a_1+a^3{x}^2+\cdots +a_{n-1}{x}^{n-2})$
所以$A(x)=A_1(x^2)+x{A}_2(x^2)$
将$x={\omega }_n^k$代入：
$$A({\omega }_n^k)=A_1({\omega }_n^{2k})+{\omega }_n^k{A}_2({\omega }_n^{2k})$$
这样分治必须保证$n=2^m$，否则无法保证分治后左右长度相同，会导致程序出错。因此我们在第一次DFT之前就要将该多项式先补成元素个数为$2^m$并且最高项次数为$n-1$的形式，补上去的项的系数为0
补成$2^m$项后，我们可以不断分治，最后每个小部分都剩下一个项，会发现一个规律。
初始: { 0 , 1 , 2 , 3 , 4 , 5 , 6 , 7 } \{0,1,2,3,4,5,6,7\} {0,1,2,3,4,5,6,7}
第一次分治： { 0 , 2 , 4 , 6 } \{0,2,4,6\} {0,2,4,6} { 1 , 3 , 5 , 7 } \{1,3,5,7\} {1,3,5,7}
第二次分治： { 0 , 4 } { 2 , 6 } { 1 , 5 } { 3 , 7 } \{0,4\}\{2,6\}\{1,5\}\{3,7\} {0,4}{2,6}{1,5}{3,7}
第三次分治： { 0 } { 4 } { 2 } { 6 } { 1 } { 5 } { 3 } { 7 } \{0\}\{4\}\{2\}\{6\}\{1\}\{5\}\{3\}\{7\} {0}{4}{2}{6}{1}{5}{3}{7}
蝶形算法
规律：
0->000->000->0,所以0在第0个位置
1->001->100->4,所以1在第4个位置
2->010->010->2,所以2在第2个位置
3->011->110->6,所以3在第6个位置
4->100->001->1,所以4在第1个位置
5->101->101->5,所以5在第5个位置
6->110->011->3,所以6在第3个位置
7->111->111->7,所以7在第7个位置
即每个数的二进制数反转后的值即为它所在位置
用该规律补全多项式后，进行傅里叶变换即可。

附kuangbin大神模板一份：

 
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <queue>
#include <cmath>
#include <string>
#include <cstring>
#include <map>
#include <set>
#include <cmath>
#include <tr1/unordered_map>
 
using namespace std;
#define me(x,y) memset(x,y,sizeof x)
#define MIN(x,y) x < y ? x : y
#define MAX(x,y) x > y ? x : y
#define re register
typedef long long ll;
typedef unsigned long long ull;
const double eps = 1e-08;
const double PI = acos(-1.0);
//复数结构体
struct Complex{
	double x,y;//实部和虚部 x+yi
	Complex(double _x = 0.0,double _y = 0.0){
		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);
	}
};
/*
* 进行 FFT 和 IFFT 前的反转变换。
* 位置 i 和(i 二进制反转后位置)互换
* len 必须为 2 的幂
*/
void change(Complex y[],int len){
	int i,j,k;
	for(i = 1, j = len/2;i <len - 1;i++){
		if(i < j)swap(y[i],y[j]);
		//交换互为小标反转的元素,i<j 保证交换一次
		//i 做正常的 +1,j 左反转类型的 +1, 始终保持 i 和 j 是反转的
		k = len/2;
		while(j >= k){
			j -= k;
			k /= 2;
		}
		if(j < k)j += k;
	}
}
/*
* 做 FFT
* len 必须为2^k形式
* on==1 时是 DFT,on==-1 时是 IDFT
*/
void fft(Complex y[],int len,int on){
	change(y,len);
	for(int h = 2; h <= len; h <<= 1){
		Complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
		for(int j = 0;j < len;j+=h){
			Complex w(1,0);
			for(int k = j;k < j+h/2;k++){
				Complex u = y[k];
				Complex t = w*y[k+h/2];
				y[k] = u+t;
				y[k+h/2] = u - t;
				w = w*wn;
			}
		}
	}
	if(on == -1)
		for(int i = 0;i < len;i++)
			y[i].x /= len;
}
const int MAXN = 200010;
Complex x1[MAXN],x2[MAXN];
char str1[MAXN/2],str2[MAXN/2];
int sum[MAXN];
int main(){
	while(scanf("%s%s",str1,str2)==2){
		int len1 = strlen(str1);
		int len2 = strlen(str2);
		int len = 1;
		while(len < len1*2 || len < len2*2)len<<=1;
		for(int i = 0;i < len1;i++)
			x1[i] = Complex(str1[len1 - 1 - i] - '0',0);
		for(int i = len1;i < len;i++)
			x1[i] = Complex(0,0);
		for(int i = 0;i < len2;i++)
			x2[i] = Complex(str2[len2 - 1 - i] - '0',0);
		for(int i = len2;i < len;i++)
			x2[i] = Complex(0,0);
		//求 DFT
		fft(x1,len,1);
		fft(x2,len,1);
		for(int i = 0;i < len;i++)
			x1[i] = x1[i]*x2[i];
		fft(x1,len, - 1);
		for(int i = 0;i < len;i++)
			sum[i] = (int)(x1[i].x+0.5);
		for(int i = 0;i < len;i++){
			sum[i+1]+=sum[i]/10;
			sum[i]%=10;
		}
		len = len1+len2 - 1;
		while(sum[len] <= 0 && len > 0)len-- ;
		for(int i = len;i >= 0;i-- )
			printf("%c",sum[i]+'0');
		printf("\n");	
	}
	return 0;
}

 
/*
 
 
*/

学fft时看过参考过的博客：
https://blog.csdn.net/qq_37136305/article/details/81184873
https://blog.csdn.net/enjoy_pascal/article/details/81478582
https://blog.csdn.net/under_sky_dxj/article/details/52778350