Skip to main content

快速傅立叶变换（FFT）

参考资料​

思路​

首先进行正变换，然后逐位相乘，最后通过逆变换得到答案。

复数​

struct Comp
{
	double real,imag;
	Comp(double real=0.0,double imag=0.0):real(real),imag(imag){}
	Comp operator+(const Comp &rhs) const{return Comp(real+rhs.real,imag+rhs.imag);}
	Comp operator-(const Comp &rhs) const{return Comp(real-rhs.real,imag-rhs.imag);}
	Comp operator*(const Comp &rhs) const{return Comp(real*rhs.real-imag*rhs.imag,real*rhs.imag+rhs.real*imag);}
	Comp operator/(const Comp &rhs) const{return Comp((real*rhs.real+imag*rhs.imag)/(rhs.real*rhs.real+rhs.imag*rhs.imag),(imag*rhs.real-real*rhs.imag)/(rhs.real*rhs.real+rhs.imag*rhs.imag));}
};

实现​

分治递归​

void fft(Comp *f,int n,int type)
{
	if(n==1)return;
	int mid=n>>1;
	Comp *g=f,*h=f+mid;
	for(int i=0;i<n;i++)t[i]=f[i];
	for(int i=0;i<mid;i++)
	{
		g[i]=t[i<<1];
		h[i]=t[i<<1|1];
	}
	fft(g,mid,type);
	fft(h,mid,type);
	Comp cur(1,0),step(cos(pi*2/n),sin(pi*2/n)*type);
	for(int i=0;i<mid;i++)
	{
		t[i]=g[i]+cur*h[i];
		t[i+mid]=g[i]-cur*h[i];
		cur=cur*step;
	}
	for(int i=0;i<n;i++)f[i]=t[i];
}

倍增迭代​

void fft(Comp *f,int n,int type)
{
	for(int i=0;i<n;i++)
	{
		if(i<r[i])swap(f[i],f[r[i]]);
	}
	for(int k=1;k<n;k<<=1)
	{
		Comp step(cos(pi/k),sin(pi/k)*type);
		for(int i=0;i<n;i+=k<<1)
		{
			Comp cur(1,0);
			for(int j=0;j<k;j++)
			{
				Comp x=f[i+j],y=f[i+j+k]*cur;
				f[i+j]=x+y;
				f[i+j+k]=x-y;
				cur=cur*step;
			}
		}
	}
}

例题​

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

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

给定一个 nn 次多项式 F(x)F(x)，和一个 mm 次多项式 G(x)G(x)。

请求出 F(x)F(x) 和 G(x)G(x) 的乘积。

Code (2)
#include <bits/stdc++.h>
using namespace std;

using Comp=complex<double>;
const double pi=acos(-1.0);
const int N=1<<20;
Comp a[N<<1],b[N<<1];
int r[N<<1];
void fft(Comp *f,int n,int type)
{
	for(int i=0;i<n;i++)if(i<r[i])swap(f[i],f[r[i]]);
	for(int k=1;k<n;k<<=1)
	{
		Comp step(cos(pi/k),sin(pi/k)*type);
		for(int i=0;i<n;i+=k<<1)
		{
			Comp cur(1,0);
			for(int j=0;j<k;j++)
			{
				Comp x=f[i+j],y=f[i+j+k]*cur;
				f[i+j]=x+y;
				f[i+j+k]=x-y;
				cur=cur*step;
			}
		}
	}
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n,m;
	cin>>n>>m;
	for(int i=0;i<=n;i++)
	{
		double x;
		cin>>x;
		a[i]=Comp(x,0);
	}
	for(int i=0;i<=m;i++)
	{
		double x;
		cin>>x;
		b[i]=Comp(x,0);
	}
	int lim=1,cnt=0;
	while(lim<=n+m){lim<<=1;cnt++;}
	for(int i=0;i<lim;i++)r[i]=r[i>>1]>>1|(i&1)<<cnt-1;
	fft(a,lim,1);
	fft(b,lim,1);
	for(int i=0;i<lim;i++)a[i]*=b[i];
	fft(a,lim,-1);
	for(int i=0;i<=n+m;i++)cout<<int(a[i].real()/lim+0.5)<<' ';
	return 0;
}

洛谷 P1919 【模板】高精度乘法 | A*B Problem 升级版​

洛谷 P1919 【模板】高精度乘法 | A*B Problem 升级版

给你两个正整数 a,ba,b，求 a×ba \times b。（1≤a,b≤1010000001\le a,b \le 10^{1000000}）

Code (2)
#include <bits/stdc++.h>
using namespace std;

using Comp=complex<double>;
const double pi=acos(-1.0);
const int N=1<<20;
Comp a[N<<1],b[N<<1];
int c[N<<1],r[N<<1];
void fft(Comp *f,int n,int type)
{
	for(int i=0;i<n;i++)if(i<r[i])swap(f[i],f[r[i]]);
	for(int k=1;k<n;k<<=1)
	{
		Comp step(cos(pi/k),sin(pi/k)*type);
		for(int i=0;i<n;i+=k<<1)
		{
			Comp cur(1,0);
			for(int j=0;j<k;j++)
			{
				Comp x=f[i+j],y=f[i+j+k]*cur;
				f[i+j]=x+y;
				f[i+j+k]=x-y;
				cur=cur*step;
			}
		}
	}
}
int main()
{
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	string s1,s2;
	cin>>s1>>s2;
	int n=s1.length()-1,m=s2.length()-1;
	for(int i=0;i<=n;i++)a[i]=Comp(s1[n-i]-'0',0);
	for(int i=0;i<=m;i++)b[i]=Comp(s2[m-i]-'0',0);
	int lim=1,cnt=0;
	while(lim<=n+m){lim<<=1;cnt++;}
	for(int i=0;i<lim;i++)r[i]=r[i>>1]>>1|(i&1)<<cnt-1;
	fft(a,lim,1);
	fft(b,lim,1);
	for(int i=0;i<lim;i++)a[i]=a[i]*b[i];
	fft(a,lim,-1);
	for(int i=0;i<=n+m;i++)c[i]=a[i].real()/lim+0.5;
	int dig=0;
	while(dig<=n+m||c[dig])
	{
		c[dig+1]+=c[dig]/10;
		c[dig++]%=10;
	}
	for(int i=dig-1;i>=0;i--)cout<<c[i];
	return 0;
}