FFT & NTT & FWT
只是学习笔记,真心推荐 cmd ,他讲的真的细到把所有的前置知识都讲了一遍。
FFT & NTT & FWT 大杂烩
首先我们引入卷积的概念,对于两个多项式进行卷积,形如:
有
其中
一般来说两个多项式进行朴素实现是
那么就下来介绍几种对于特殊的卷积可以加快速度的算法。
一、FFT
FFT 是针对于加法卷积(也就是多项式相乘)的快速计算方法。
大体思路
首先我们知道对于一个
但是朴素求值和插值都是
加速
前置:单位根
首先大家应该学过复数,接下来介绍单位根 。
有一个单位圆,我们通过
如图是
有一天FFT的发明者 傅里叶 突然想把
首先由几个性质在几何角度看来非常容易证明:
1、
2、
3、
4、
5、
6、
主体
那么接下来讲解 FFT 的主要过程:
求值:
对一个长度为
对于
那么我们尝试写出
接下来就是单位根出场的时候了,我们将
对于
对于
上面的式子意味着什么,它意味着我们可以通过分治去
好,我们现在已经可以
插值:
设我们刚才把
算了还是证明一下:
我们将
分类讨论,对于
贡献为
对于
贡献为
所以最终就有
所以我们将插值也转成了求值,同样可以
实现
朴素实现
现在可以再来看一开始的流程图
整个过程就比较清晰了吧。
给出模板题: P3803 【模板】多项式乘法(FFT)
码(未卡常版本)
#include<bits/stdc++.h> using namespace std; typedef long long ll; typedef pair<int,int> pii; typedef unsigned long long ull; #define mk make_pair #define ps push_back #define fi first #define se second const int N=5e6+10,inf=0x3f3f3f3f; const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f; inline ll read(){ char c=getchar();ll x=0;bool f=0; while(!isdigit(c))f=c=='-'?1:0,c=getchar(); while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar(); return f?-x:x; } const double PI=acos(-1); int n,m; struct jj{ // 复数结构体 double x,y; inline jj operator +(const jj&k){return {x+k.x,y+k.y};} inline jj operator -(const jj&k){return {x-k.x,y-k.y};} inline jj operator *(const jj&k){return {x*k.x-y*k.y,x*k.y+y*k.x};} }f[N],g[N],tp[N]; inline void fft(jj *f,int l,int r,bool fl){ if(l==r)return; int mid=l+r>>1,len=r-l+1; for(int i=l;i<=r;++i) tp[i]=f[i]; for(int i=l;i<=mid;++i){ f[i]=tp[(i-l)*2+l],f[mid+1+(i-l)]=tp[(i-l)*2+l+1];// 按照奇偶分裂 } fft(f,l,mid,fl),fft(f,mid+1,r,fl);// 分治计算 jj op={cos(2*PI/len),sin(2*PI/len)},now={1,0};// op=w_n^1 if(fl)op.y*=-1;// 如果 fl 为 1 表示是插值,此时 op=w_n^-1 for(int i=l;i<=mid;++i){ tp[i]=f[i]+now*f[mid+1+(i-l)]; tp[mid+1+(i-l)]=f[i]-now*f[mid+1+(i-l)]; now=now*op; } for(int i=l;i<=r;++i) f[i]=tp[i]; } signed main(){ #ifndef ONLINE_JUDGE freopen("in.in","r",stdin); freopen("out.out","w",stdout); #endif ios::sync_with_stdio(0),cin.tie(0),cout.tie(0); n=read(),m=read(); for(int i=0;i<=n;++i) f[i]={read(),0}; for(int i=0;i<=m;++i) g[i]={read(),0}; m+=n;n=1; while(n<=m)n<<=1; --n; fft(f,0,n,0),fft(g,0,n,0); for(int i=0;i<=n;++i) f[i]=f[i]*g[i]; fft(f,0,n,1); for(int i=0;i<=m;++i){ cout<<(int)(f[i].x/(n+1)+0.499)<<' ';//有精度问题,四舍五入 } }
常数非常大,luogu最大点跑了 800 ms ,考虑一些优化。
优化
首先可以把递归改为非递归版的迭代。其次一个重要优化是 蝴蝶变换 。他说的是我们递归的过程会把
首先证明这个事:
我们脑动模拟一下递归过程,每次按照偶在左半边,奇在右半边去分,相当于看二进制下最后一位,按照大小确定了最终二进制下第一位的大小。
有点抽象,但是有图:
以下将经过
这是一开始的
然后我们按照最后一位的大小,分开
这时可以看到
那么下一层也就根据
而事实也确实如此。
好,设 pos[i]=(pos[i>>1]>>1)|(i&1?n>>1:0)
,
所以我们就获得了一个常数较小的 FFT:
码(卡常版)
#include<bits/stdc++.h> using namespace std; typedef long long ll; typedef pair<int,int> pii; typedef unsigned long long ull; #define mk make_pair #define ps push_back #define fi first #define se second const int N=1e6+10,inf=0x3f3f3f3f; const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f; inline ll read(){ char c=getchar();ll x=0;bool f=0; while(!isdigit(c))f=c=='-'?1:0,c=getchar(); while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar(); return f?-x:x; } const double PI=acos(-1); struct jj{ double x,y; inline jj operator +(const jj&k){return {x+k.x,y+k.y};} inline jj operator -(const jj&k){return {x-k.x,y-k.y};} inline jj operator *(const jj&k){return {x*k.x-y*k.y,x*k.y+y*k.x};} }f[N<<2],g[N<<2],ji[N<<2]; int pos[N<<2]; inline void fft(jj *f,int n,bool fl){ for(int i=0;i<n;++i) if(i<pos[i])swap(f[i],f[pos[i]]);//pos[i] 表示预处理二进制翻转后的位置 for(int len=2,mid;len<=n;len<<=1){ mid=len>>1; jj op=ji[len],lp;//ji[len]表示预处理 w_n^1 ,因为 sin,cos 太慢了,所以预处理会少调用几次 if(fl)op.y*=-1; for(int i=0;i<n;i+=len){ jj now={1,0}; for(int j=i;j<i+mid;++j){ lp=now*f[j+mid]; f[j+mid]=f[j]-lp;f[j]=f[j]+lp; now=now*op; } } } } int n,m; signed main(){ #ifndef ONLINE_JUDGE freopen("in.in","r",stdin); freopen("out.out","w",stdout); #endif ios::sync_with_stdio(0),cin.tie(0),cout.tie(0); n=read(),m=read(); for(int i=0;i<=n;++i) f[i]={read(),0}; for(int i=0;i<=m;++i) g[i]={read(),0}; m+=n;n=1;ji[1]={1,0}; while(n<=m)n<<=1,ji[n]={cos(2*PI/n),sin(2*PI/n)}; for(int i=0;i<n;++i) pos[i]=(pos[i>>1]>>1)|(i&1?n>>1:0); fft(f,n,0),fft(g,n,0); for(int i=0;i<=n;++i) f[i]=f[i]*g[i]; fft(f,n,1); for(int i=0;i<=m;++i){ cout<<(int)(f[i].x/n+0.499)<<' '; } }
luogu最大点 469 ms ,而且出奇的好写。
upd:经 Qyun 和 CuFeO4 推荐更新“三次变两次优化”
三次变两次
因为是复数运算,设一个复数多项式
所以我们要求的就是
码(最终版)
#include<bits/stdc++.h> using namespace std; typedef long long ll; typedef pair<int,int> pii; typedef unsigned long long ull; #define mk make_pair #define ps push_back #define fi first #define se second const int N=1e6+10,inf=0x3f3f3f3f; const ll mod=1e9+7,linf=0x3f3f3f3f3f3f3f3f; inline ll read(){ char c=getchar();ll x=0;bool f=0; while(!isdigit(c))f=c=='-'?1:0,c=getchar(); while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar(); return f?-x:x; } const double PI=acos(-1); struct jj{ double x,y; inline jj operator +(const jj&k){return {x+k.x,y+k.y};} inline jj operator -(const jj&k){return {x-k.x,y-k.y};} inline jj operator *(const jj&k){return {x*k.x-y*k.y,x*k.y+y*k.x};} }f[N<<2],ji[N<<2]; int pos[N<<2]; inline void fft(jj *f,int n,bool fl){ for(int i=0;i<n;++i) if(i<pos[i])swap(f[i],f[pos[i]]); for(int len=2,mid;len<=n;len<<=1){ mid=len>>1; jj op=ji[len],lp; if(fl)op.y*=-1; for(int i=0;i<n;i+=len){ jj now={1,0}; for(int j=i;j<i+mid;++j){ lp=now*f[j+mid]; f[j+mid]=f[j]-lp;f[j]=f[j]+lp; now=now*op; } } } } int n,m; signed main(){ #ifndef ONLINE_JUDGE freopen("in.in","r",stdin); freopen("out.out","w",stdout); #endif ios::sync_with_stdio(0),cin.tie(0),cout.tie(0); n=read(),m=read(); for(int i=0;i<=n;++i) f[i].x=read(); for(int i=0;i<=m;++i) f[i].y=read(); m+=n;n=1;ji[1]={1,0}; while(n<=m)n<<=1,ji[n]={cos(2*PI/n),sin(2*PI/n)}; for(int i=0;i<n;++i) pos[i]=(pos[i>>1]>>1)|(i&1?n>>1:0); fft(f,n,0); for(int i=0;i<=n;++i) f[i]=f[i]*f[i]; fft(f,n,1); for(int i=0;i<=m;++i){ cout<<(int)(f[i].y/n/2+0.499)<<' '; } }
luogu最大点 348 ms。
例题
P1919 【模板】高精度乘法 | A*B Problem 升级版
另一个板子题,FFT 加速即可。
P3338 [ZJOI2014] 力
给出序列
求所有的
我们考虑如何把
前面已经是卷积的形式了,无需转化,后面的不是,我们考虑给他转成卷积的形式。
画个图看看:
他们是这样相乘的,但是卷积是这样相乘的:
所以我们考虑给
对两者 FFT 然后加起来即可。
二、 NTT
NTT ,中文 快速数论变换 是用来解决卷积过程中需要取模的快速加法卷积,说白了就是可以取模的 FFT。
这个东西就是直接把
码
#include<bits/stdc++.h> using namespace std; typedef long long ll; #define int ll typedef pair<int,int> pii; typedef unsigned long long ull; #define mk make_pair #define ps push_back #define fi first #define se second const int N=1e6+10,inf=0x3f3f3f3f; const ll mod=998244353,linf=0x3f3f3f3f3f3f3f3f,g=3; inline ll read(){ char c=getchar();ll x=0;bool f=0; while(!isdigit(c))f=c=='-'?1:0,c=getchar(); while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar(); return f?-x:x; } inline ll qpow(ll x,ll y){ ll ans=1; while(y){ if(y&1)ans=ans*x%mod; x=x*x%mod;y>>=1; } return ans; } ll ng=qpow(g,mod-2);// g 的逆元 int pos[N<<2]; inline void ntt(ll *f,int n,bool fl){ for(int i=0;i<n;++i) if(i<pos[i])swap(f[i],f[pos[i]]); for(int len=2,mid;len<=n;len<<=1){ int op=qpow(fl?g:ng,(mod-1)/len),tp;mid=len>>1; for(int i=0;i<n;i+=len){ ll now=1; for(int j=i;j<i+mid;++j){ tp=f[j+mid]*now%mod; f[j+mid]=f[j]-tp; if(f[j+mid]<0)f[j+mid]+=mod; f[j]+=tp; if(f[j]>=mod)f[j]-=mod; now=now*op%mod; } } } } int n,m; ll f[N<<2],ff[N<<2]; signed main(){ #ifndef ONLINE_JUDGE freopen("in.in","r",stdin); freopen("out.out","w",stdout); #endif ios::sync_with_stdio(0),cin.tie(0),cout.tie(0); n=read(),m=read(); for(int i=0;i<=n;++i) f[i]=read(); for(int i=0;i<=m;++i) ff[i]=read(); m+=n;n=1; while(n<=m)n<<=1; for(int i=0;i<n;++i) pos[i]=(pos[i>>1]>>1)|(i&1?n>>1:0); ntt(f,n,0),ntt(ff,n,0); for(int i=0;i<n;++i) f[i]=f[i]*ff[i]%mod; ntt(f,n,1); ll ny=qpow(n,mod-2); for(int i=0;i<=m;++i) cout<<f[i]*ny%mod<<' '; }
但是我的 NTT 好像没有 FFT 快???
三、 FWT
在 OI 中,FWT是用来解决有关位运算卷积的快速卷积算法。
即对于
或运算
此时
我们设
考虑如何快速求
注意到将
再看
其中
类似的可以退出逆运算为
也就可以在得到
注
感觉就和 FFT , NTT 一样, FWT(|) 找了一种特殊运算,使得可以分治
与运算
类比或运算,这里直接给出式子
求值:
插值:
异或运算
异或运算是 FWT 中最难理解的一个,是一种非常神奇的构造。
以下
首先我们定义
证明:
我们一位一位的去考虑,在一开始,两边都是 0,然后考虑新的一位,如果要做出改变,首先
综上,
接下来我们定义
接下来考虑如何快速计算
首先举一个例子:
假设
然后我们一分为二,递归求解,但是注意递归到下一层后,他们在本层的第一位将不再考虑,也就是
这也就是为什么之前注里面说分治的过程是逐渐考虑二进制位的影响。
理解了这里,下面就很好理解了。
对于
同时有逆推
模板题
P4717 【模板】快速莫比乌斯/沃尔什变换 (FMT/FWT)
码
#include<bits/stdc++.h> using namespace std; typedef long long ll; #define int ll typedef pair<int,int> pii; typedef unsigned long long ull; #define mk make_pair #define ps push_back #define fi first #define se second const int N=1e6+10,inf=0x3f3f3f3f; const ll linf=0x3f3f3f3f3f3f3f3f,mod=998244353; inline ll read(){ char c=getchar();ll x=0;bool f=0; while(!isdigit(c))f=c=='-'?1:0,c=getchar(); while(isdigit(c))x=(x<<1)+(x<<3)+(c^48),c=getchar(); return f?-x:x; } inline void OR(int *f,int n,int op){ for(int len=2;len<=n;len<<=1){ int mid=len>>1; for(int j=0;j<n;j+=len){ for(int k=j;k<j+mid;++k){ f[k+mid]=(f[k+mid]+f[k]*op)%mod; } } } } inline void AND(int *f,int n,int op){ for(int len=2;len<=n;len<<=1){ int mid=len>>1; for(int j=0;j<n;j+=len){ for(int k=j;k<j+mid;++k) f[k]=(f[k]+f[k+mid]*op)%mod; } } } inline void XOR(int *f,int n,int op){ for(int len=2;len<=n;len<<=1){ int mid=len>>1,tp; for(int j=0;j<n;j+=len){ for(int k=j;k<j+mid;++k){ tp=f[k]; f[k]=(f[k]+f[k+mid])*op%mod; f[k+mid]=(tp-f[k+mid])*op%mod; } } } } inline ll qpow(ll x,ll y){ ll ans=1; while(y){ if(y&1)ans=ans*x%mod; x=x*x%mod;y>>=1; } return ans; } int a[N],b[N],c[N],n; signed main(){ #ifndef ONLINE_JUDGE freopen("in.in","r",stdin); freopen("out.out","w",stdout); #endif ios::sync_with_stdio(0),cin.tie(0),cout.tie(0); n=read();int m=1<<n; for(int i=0;i<m;++i) a[i]=read(); for(int i=0;i<m;++i) b[i]=read(); OR(a,m,1),OR(b,m,1); for(int i=0;i<m;++i) c[i]=a[i]*b[i]%mod; OR(c,m,-1);OR(a,m,-1),OR(b,m,-1); for(int i=0;i<m;++i) cout<<(c[i]+mod)%mod<<' '; cout<<'\n'; AND(a,m,1);AND(b,m,1); for(int i=0;i<m;++i) c[i]=a[i]*b[i]%mod; AND(a,m,-1),AND(b,m,-1),AND(c,m,-1); for(int i=0;i<m;++i) cout<<(c[i]+mod)%mod<<' '; cout<<'\n'; XOR(a,m,1),XOR(b,m,1); for(int i=0;i<m;++i) c[i]=a[i]*b[i]%mod; XOR(c,m,qpow(2,mod-2)); for(int i=0;i<m;++i) cout<<(c[i]+mod)%mod<<' '; }
小结
感觉三者都是通过构造找了一些特殊的点或特殊的运算,使得我们可以通过这写特殊的性质达到快速计算的效果,当然这一切的前提都是因为我们用到的
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 2分钟学会 DeepSeek API,竟然比官方更好用!
· .NET 使用 DeepSeek R1 开发智能 AI 客户端
· 10亿数据,如何做迁移?
· 推荐几款开源且免费的 .NET MAUI 组件库
· c# 半导体/led行业 晶圆片WaferMap实现 map图实现入门篇