快速傅里叶变换学习笔记

一.单位根

满足 $x^n-1=0$ 的 $x$ 称为 $n$ 次单位根

1.本原单位根

若 $\omega$ 是 $n$ 次单位根,且 $\omega^0$ , $\omega^1$ , $\omega^2$ ,..., $\omega^{n-1}$ 恰好是所有的 $n$ 次单位根，则称 $\omega$ 为 $n$ 次本原单位根,记作 $\omega_n$

2.单位根的计算

在复数域 $\mathbb{C}$ 上， $e^{ix}=cosx+isinx$ （欧拉公式，证明的话左边将 $\exp$ (ix) 展开成 $\texttt{EGF}$ 的形式，右边同理即证）

重要的性质: $n$ 次单位根在复平面上等分单位圆

由"重要的性质"可以得到, $\omega_n$ = $\exp(\frac{2 \pi i}{n})$ = $cos \frac{2\pi }{n}$ + $isin \frac{2\pi }{n}$

以此我们可以表示出所有的 $n$ 次单位根，即: $\omega_n^k$ = $cos \frac{2\pi k}{n}$ + $isin \frac{2\pi k}{n}$

二.离散傅里叶变换DFT

前置芝士1:多项式的点值表达

对于一个多项式 $f(x)$ ，我们代入每个 $x_i$ ,会得到对应的 $f(x_i)$ ,这些 $(x_i,f(x_i))$ 构成了这个多项式的点值表达，且唯一确定了这个多项式

前置芝士2：多项式相乘

两个多项式的乘积被定义为:
$A(x)B(x)$ = $\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{n}{a_ib_jx^{i+j}}$ = $\sum\limits_{i=1}^{2n}{c_kx^k}$
其中 $c$ 是 $a$ 和 $b$ 的卷积,计算两个多项式相乘的朴素算法是 $O(n^2)$ 的
如果将两个多项式先转换为点值表达,再逐项相乘，时间复杂度只有 $O(n)$
即 $(x,A(x))*(x,B(x))=(x,A(x)*B(x)）$
那么，复杂度优化的瓶颈在于将多项式转化为点值表达（以及转回去），FFT解决的正是这个问题

1.未经优化的DFT

设 $a$ 是长度为 $n$ 的数列，对于 $0\leq k <n$ ,定义

$b_k$ = $\sum\limits_{i=0}^{n-1}{(a_i\omega^{ki})}$

其中 $b$ 即为 $a$ 的 $DFT$

这是，我们令多项式 $f(x)$ = $\sum a_ix^i$ ,则 $b_k$ 就是 $f(x)$ 在 $\omega^k$ 处的点值

当我们计算完两个多项式乘积的点值,如何将其转回多项式?

2.DFT的逆变换IDFT

$a_k$ = $\frac{1}{n}\sum\limits_{i=0}^{n-1}{b_i\omega^{-ki}}$

这个柿子与DFT的相似度极高的特性使得我们不需要重新写一个函数来处理IDFT，只需要提前预处理单位根的倒数（即共轭复数），利用FFT转化，最后全部除以 $n$ 即可

四.快速傅里叶变换FFT

前置芝士:单位根的性质

$(1)\omega_{2n}^{2k}=\omega_n^k$

$(2)\omega_{2n}^{n+k}=- \omega_{2n}^k$

1.FFT

啥是FFT？就是利用DFT的奇偶性质快速计算DFT的一个分治算法

设 $n=2m$ ,将 $A(x)$ 按次数奇偶分类，有:

$A(x)$ = $\sum\limits_{i=0}^{n-1}a_ix^i$

= $\sum\limits_{i=0}^{m-1}a_{2i}x^{2i}$ + $\sum\limits_{i=0}^{m-1}a_{2i+1}x^{2i+1}$

= $\sum\limits_{i=0}^{m-1}a_{2i}x^{2i}$ + $x\sum\limits_{i=0}^{m-1}a_{2i+1}x^{2i}$

我们令 $A_0(x)=\sum\limits_{i=0}^{m-1}a_{2i}x^{i}$ , $A_1(x)=\sum\limits_{i=0}^{m-1}a_{2i+1}x^{i}$

则 $A(x)=A_0(x^2)+xA_1(x^2)$

通过上式我们可以得出，如果我们得到了 $A_0(x)$ $A_1(X)$ 在 $\omega_m^0$ , $\omega_m^1$ , $\omega_m^2$ ,..., $\omega_m^{n-1}$ 处的点值，则可在 $O(n)$ 内得到 $A(x)$ 在 $\omega_m^0$ , $\omega_m^1$ , $\omega_m^2$ ,..., $\omega_m^{n-1}$ 处的点值

而 $A_0(x),A_1(x)$ 是可以分治计算的,递归深度不会超过 $\log n$

综上，我们可以在 $O(n \log n)$ 内完成对 $A(x)$ 的点值求值

五.一些优化

1.位逆序置换

我们发现，令 $rev(x)$ 表示 $x$ 经过二进制反转后的数，且令 $b_i=a_{rev(i)}$ ,则每次对 $a_i$ 进行的操作在 $b_i$ 中变为了对相邻两个序列进行的操作，那么我们只需要预先处理出 $b_i$ ，直接递归向上不断还原即可
它在代码中一般这样子体现：

fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));

2.蝴蝶操作

因为double的乘法比较慢，我们在代码中有这样一段：

Complex x=arr[j+k];
arr[j+k]=x+w*arr[j+mid+k],arr[j+mid+k]=x-w*arr[j+mid+k];

我们发现w*arr[j+mid+k]这东西被算了两次，所以应该这么写

Complex x=arr[j+k],y=w*arr[j+mid+k];
arr[j+k]=x+y,arr[j+mid+k]=x-y;

然后就优化完了/cy
~~所以这难道不是常数优化？~~

六.代码实现

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

#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
	typedef long long ll;
	typedef double db;
	#define pf printf
	#define pc putchar
	#define fr(i,x,y) for(register int i=(x);i<=(y);++i)
	#define pfr(i,x,y) for(register int i=(x);i>=(y);--i)
	#define go(u) for(int i=head[u];i;i=e[i].nxt)
	#define enter pc('\n')
	#define space pc(' ')
	#define fir first
	#define sec second
	#define MP make_pair
	const int inf=0x3f3f3f3f;
	const ll inff=1e15;
	inline int read()
	{
		int sum=0,f=1;
		char ch=0;
		while(!isdigit(ch))
		{
			if(ch=='-') f=-1;
			ch=getchar();
		}
		while(isdigit(ch))
		{
			sum=sum*10+(ch^48);
			ch=getchar();
		}
		return sum*f;
	}
	inline void write(int x)
	{
		if(x<0)
		{
			x=-x;
			pc('-');
		}
		if(x>9) write(x/10);
		pc(x%10+'0');
	}
	inline void writeln(int x)
	{
		write(x);
		enter;
	}
	inline void writesp(int x)
	{
		write(x);
		space;
	}
}
using namespace my_std;
const int maxn=1e7+50;
struct Complex
{
	double x,y;
	Complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
double PI=acos(-1.0);
Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int N,M,limit=1,l,r[maxn];
inline void fft(Complex *arr,int pd)
{
	fr(i,0,limit-1) if(i<r[i]) swap(arr[i],arr[r[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		Complex Wn(cos(PI/mid),pd*sin(PI/mid));
		for(int j=0,R=mid<<1;j<limit;j+=R)
		{
			Complex w(1,0);
			for(int k=0;k<mid;++k,w=w*Wn)
			{
				Complex x=arr[j+k],y=w*arr[j+mid+k];
				arr[j+k]=x+y,arr[j+mid+k]=x-y; 
			} 
		}
	}
}
int main(void)
{
	N=read(),M=read();
	fr(i,0,N) a[i].x=read();
	fr(i,0,M) b[i].x=read();
	//system("pause");
	while(limit<=M+N) limit<<=1,++l;
	fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	fft(a,1),fft(b,1);
	fr(i,0,limit) a[i]=a[i]*b[i];
	fft(a,-1);
	fr(i,0,M+N) writesp((int)(a[i].x/limit+0.5));
	return 0;
}

快速数论变换NTT

#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
	typedef long long ll;
	typedef double db;
	#define pf printf
	#define pc putchar
	#define fr(i,x,y) for(register ll i=(x);i<=(y);++i)
	#define pfr(i,x,y) for(register ll i=(x);i>=(y);--i)
	#define go(u) for(ll i=head[u];i;i=e[i].nxt)
	#define enter pc('\n')
	#define space pc(' ')
	#define fir first
	#define sec second
	#define MP make_pair
	const ll inf=0x3f3f3f3f;
	const ll inff=1e15;
	inline ll read()
	{
		ll sum=0,f=1;
		char ch=0;
		while(!isdigit(ch))
		{
			if(ch=='-') f=-1;
			ch=getchar();
		}
		while(isdigit(ch))
		{
			sum=sum*10+(ch^48);
			ch=getchar();
		}
		return sum*f;
	}
	inline void write(ll x)
	{
		if(x<0)
		{
			x=-x;
			pc('-');
		}
		if(x>9) write(x/10);
		pc(x%10+'0');
	}
	inline void writeln(ll x)
	{
		write(x);
		enter;
	}
	inline void writesp(ll x)
	{
		write(x);
		space;
	}
}
using namespace my_std;
const ll maxn=1e7+50,G=3,mod=998244353;
ll N,M,limit=1,a[maxn],b[maxn],l,r[maxn];
inline ll ksmod(ll a,ll b)
{
	ll ans=1;
	while(b)
	{
		if(b&1) ans=(ans*a)%mod;
		a=(a*a)%mod;
		b>>=1;
	}
	return ans%mod;
}
inline void NTT(ll *a,ll pd)
{
	fr(i,0,limit-1) if(i<r[i]) swap(a[i],a[r[i]]);
	for(ll i=1;i<limit;i<<=1)
	{
		ll gn=ksmod(G,(mod-1)/(i<<1));
		for(ll j=0;j<limit;j+=(i<<1))
		{
			ll g=1;
			for(ll k=0;k<i;++k,g=(g*gn)%mod)
			{
			 	ll x=a[j+k],y=g*a[j+k+i]%mod;
				a[j+k]=(x+y)%mod,a[j+k+i]=(x-y+mod)%mod;
			}
		}
	}
	if(pd==1) return ;
	ll inv=ksmod(limit,mod-2);
	reverse(a+1,a+limit);
	fr(i,0,limit-1) a[i]=a[i]*inv%mod;
}
int main(void)
{
	N=read(),M=read();
	fr(i,0,N) a[i]=read();
	fr(i,0,M) b[i]=read();
	while(limit<=M+N) limit<<=1,++l;
	fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	NTT(a,1),NTT(b,1);
	fr(i,0,limit-1) a[i]=(a[i]*b[i])%mod;
	NTT(a,-1);
	fr(i,0,M+N) writesp(a[i]%mod);
	return 0;
}

洛谷1919【模板】A*B Problem升级版（FFT快速傅里叶）

#include<bits/stdc++.h>
using namespace std;
namespace my_std
{
	typedef long long ll;
	typedef double db;
	#define pf printf
	#define pc putchar
	#define fr(i,x,y) for(register int i=(x);i<=(y);++i)
	#define pfr(i,x,y) for(register int i=(x);i>=(y);--i)
	#define go(u) for(int i=head[u];i;i=e[i].nxt)
	#define enter pc('\n')
	#define space pc(' ')
	#define fir first
	#define sec second
	#define MP make_pair
	const int inf=0x3f3f3f3f;
	const ll inff=1e15;
	inline int read()
	{
		int sum=0,f=1;
		char ch=0;
		while(!isdigit(ch))
		{
			if(ch=='-') f=-1;
			ch=getchar();
		}
		while(isdigit(ch))
		{
			sum=sum*10+(ch^48);
			ch=getchar();
		}
		return sum*f;
	}
	inline void write(int x)
	{
		if(x<0)
		{
			x=-x;
			pc('-');
		}
		if(x>9) write(x/10);
		pc(x%10+'0');
	}
	inline void writeln(int x)
	{
		write(x);
		enter;
	}
	inline void writesp(int x)
	{
		write(x);
		space;
	}
}
using namespace my_std;
const int maxn=2e6+50;
struct Complex
{
	double x,y;
	Complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[maxn],b[maxn];
double PI=acos(-1.0);
Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int N,M,limit=1,l,r[maxn],ans[maxn];
char sa[maxn],sb[maxn]; 
inline void fft(Complex *arr,int pd)
{
	fr(i,0,limit-1) if(i<r[i]) swap(arr[i],arr[r[i]]);
	for(int mid=1;mid<limit;mid<<=1)
	{
		Complex Wn(cos(PI/mid),pd*sin(PI/mid));
		for(int j=0,R=mid<<1;j<limit;j+=R)
		{
			Complex w(1,0);
			for(int k=0;k<mid;++k,w=w*Wn)
			{
				Complex x=arr[j+k],y=w*arr[j+mid+k];
				arr[j+k]=x+y,arr[j+mid+k]=x-y; 
			} 
		}
	}
}
int main(void)
{
	scanf("%s%s",sa,sb);
	int lena=0,lenb=0,hhh=strlen(sa),hhhh=strlen(sb);
	pfr(i,hhh-1,0) a[lena++].x=sa[i]-48;
	pfr(i,hhhh-1,0) b[lenb++].x=sb[i]-48;
	while(limit<lena+lenb) limit<<=1,++l;
	fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	fft(a,1),fft(b,1);
	fr(i,0,limit) a[i]=a[i]*b[i];
	fft(a,-1);
	int tot=0;
	fr(i,0,limit)
	{
		ans[i]+=(int)(a[i].x/limit+0.5);
		if(ans[i]>=10) ans[i+1]+=ans[i]/10,ans[i]%=10,limit+=(i==limit);
	}
	while(!ans[limit]&&limit>=1) --limit;
	++limit;
	while(--limit>=0) write(ans[limit]);
	return 0; 
}

完结撒花！！！

上一篇多项式牛顿迭代学习笔记

下一篇2020提高组十连测day2 毒瘤的区间题解

本文作者：L_G_J

本文链接：https://www.cnblogs.com/lgj-lgj/p/12230262.html

posted @ 2020-07-24 07:56 L_G_J 阅读(414) 评论(3) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 分享一个免费、快速、无限量使用的满血 DeepSeek R1 模型，支持深度思考和联网搜索！
· 基于 Docker 搭建 FRP 内网穿透开源项目（很简单哒）
· ollama系列1：轻松3步本地部署deepseek，普通电脑可用
· 按钮权限的设计及实现
· 25岁的心里话

蒟蒻L_G_J的blog

逃げてばかりは勝てない、絶対！

快速傅里叶变换学习笔记

一.单位根

1.本原单位根

2.单位根的计算

重要的性质: $n$ 次单位根在复平面上等分单位圆

二.离散傅里叶变换DFT

前置芝士1:多项式的点值表达

前置芝士2：多项式相乘

1.未经优化的DFT

2.DFT的逆变换IDFT

四.快速傅里叶变换FFT

前置芝士:单位根的性质

1.FFT

五.一些优化

1.位逆序置换

2.蝴蝶操作

六.代码实现

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

快速数论变换NTT

洛谷1919【模板】A*B Problem升级版（FFT快速傅里叶）

公告

常用链接

最新随笔

我的标签

随笔档案 (54)

阅读排行榜

推荐排行榜

蒟蒻L_G_J的blog

逃げてばかりは勝てない、絶対！

快速傅里叶变换学习笔记

一.单位根

1.本原单位根

2.单位根的计算

重要的性质:nn次单位根在复平面上等分单位圆

二.离散傅里叶变换DFT

前置芝士1:多项式的点值表达

前置芝士2：多项式相乘

1.未经优化的DFT

2.DFT的逆变换IDFT

四.快速傅里叶变换FFT

前置芝士:单位根的性质

1.FFT

五.一些优化

1.位逆序置换

2.蝴蝶操作

六.代码实现

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

快速数论变换NTT

洛谷1919【模板】A*B Problem升级版（FFT快速傅里叶）

公告

常用链接

最新随笔

我的标签

随笔档案 (54)

阅读排行榜

推荐排行榜

重要的性质: $n$ 次单位根在复平面上等分单位圆