FFT & NTT & FWT

只是学习笔记,真心推荐 cmd ,他讲的真的细到把所有的前置知识都讲了一遍。

\[FBI \ WARNING:本篇 NTT 部分非常不完善 \]

FFT & NTT & FWT 大杂烩

首先我们引入卷积的概念,对于两个多项式进行卷积,形如:

\[F(x) = f(x) \oplus g(x) \]

\[F_n = \sum_{i \oplus j = n} f_i g_j \]

其中 $ F_n , f_i , g_j $ 均表示多项式 $ F , f, g $ 的某一项系数。

一般来说两个多项式进行朴素实现是 $ O(mn) $ , $ m ,n $ 分别是两个多项式的次数。

那么就下来介绍几种对于特殊的卷积可以加快速度的算法。


一、FFT

FFT 是针对于加法卷积(也就是多项式相乘)的快速计算方法。

大体思路

首先我们知道对于一个 $ n $ 次多项式可以由 $ n+1 $ 个不同的点来表示,而对于两个多项式 $ f(x) , g(x) $ 相乘,我们可以通过 $ n+1 $ 个点在 $ O(n) $ 时间复杂度内完成,然后在通过插值将得到的 $ n+1 $ 个点还原成一个多项式。这就是 FFT 的大体思路。

image

但是朴素求值和插值都是 $ O(n^2) $ 的,所以我们现在要优化求值和插值的过程。

加速

前置:单位根

首先大家应该学过复数,接下来介绍单位根

有一个单位圆,我们通过 $ n $ 条线将它均分成了 $ n $ 份,我们依次编号 $ 0,1,2,\cdots ,n-1 $ ,他们对应的复数就是单位根 $ w_n^{i} $

image

如图是 $ n = 4 $ 的情况,他恰好对应着坐标轴。

有一天FFT的发明者 傅里叶 突然想把 $ w_n^i $ 代入到多项式中,然后发现这玩意有很好的性质以至于可以快速求值和插值。

首先由几个性质在几何角度看来非常容易证明:

1、 $ w_n^0 = 1 $

2、 $ w_n^k = w_n^{ k \mod n } $

3、 $ (w_n^k) ^j = w_n^{jk} $

4、 $ w_n^k \times w_n^j = w_n^{k+j} $

5、 $ w_{2n}^{2k} = w_n^k $

6、 $ n $ 是偶数, $ w_n^{k + \frac{n}{2}} = - w_n^k $


主体

那么接下来讲解 FFT 的主要过程:

求值:

对一个长度为 $ n $ 的序列,它对应着一个 $ n-1 $ 次多项式 $ F(x) $ 。

对于 $ F(x) $ ,我们按奇偶把他一分为二,下标为偶数的在前面,称为 $ Fl(x) $ ,下标为奇数的在后面,称为 $ Fr(x) $ 。

image

那么我们尝试写出 $ F(x) $ 与 $ Fl(x) , Fr(x) $ 的关系式:

\[F(x) = Fl(x^2) + x Fr(x^2) \]

接下来就是单位根出场的时候了,我们将 $ w_n^k $ 代入 $ F(x) $ 和上述式子:

对于 $ k \lt \frac{n}{2} $ 有

\[\begin{aligned} F( w_n^k ) &= Fl( w_n^{2k} ) + w_n^k Fr( w_n^{2k} ) \\ &= Fl(w_{\frac{n}{2}}^k) + w_n^k Fr( w_{\frac{n}{2}}^k ) \end{aligned} \]

对于 $ k \lt \frac{n}{2} $ ,我们代入 $ w_n^{k + \frac{n}{2}} $ 有

\[\begin{aligned} F(w_n^{k+\frac{n}{2}}) &= Fl(w_n^{2k+n}) + w_n^{k+\frac{n}{2}} Fr(w_n^{2k+n}) \\ &= Fl(w_{\frac{n}{2}}^k) - w_n^k Fr(w_{\frac{n}{2}}^k) \end{aligned} \]

上面的式子意味着什么,它意味着我们可以通过分治去 $ O(n \log n) $ 计算 $ F(w_n^k) $ 的值。

好,我们现在已经可以 $ O(n \log n) $ 的将一个多项式转化为 $ ( w_n^k , F( w_n^k ) ) $ 的点值表达式了,然后可以 $ O(n) $ 点乘出最终的多项式 , 然后我们现在需要快速将这个多项式插回去。

插值:

设我们刚才把 $ F(x) $ 转成了 $ G(x) $ 。那么我们直接不加证明的给出:

\[n F(k) = \sum_{i=0}^{n-1} (w_n^{-k})^i G(i) \]

算了还是证明一下:

我们将 $ G(i) $ 代入

\[\begin{aligned} 右边 &= \sum_{i=0}^{n-1} w_n^{-ik} \sum_{j=0}^{n-1} w_n^{ij} F(j) \\ &= \sum_{j=0}^{n-1} F(j) w_n^{ij} \sum_{i=0}^{n-1} w_n^{-ik} \\ &= \sum_{j=0}^{n-1} F(j) \sum_{i=0}^{n-1} w_n^{i(j-k)} \end{aligned} \]

分类讨论,对于 $ j=k $ 的情况:

贡献为 $ F(k) \sum_{i=0}^{n-1} w_n^0 = n F(k) $

对于 $ j \not = k $ 的情况:

贡献为 $ F(j) \sum_{i=0}^{n-1} w_n^{i p} $ , 等比数列求和, $ \sum_{i=0}^{n-1} w_n^{i p} = \frac{ w_n^{ n p } - 1 }{ w_n^p - 1 } = 0 $ ,也就是没有贡献!

所以最终就有

\[右边 = n F(k) = 左边 \]

所以我们将插值也转成了求值,同样可以 $ O(n \log n) $ 求。

实现

朴素实现

现在可以再来看一开始的流程图

image

整个过程就比较清晰了吧。

给出模板题: 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 ,考虑一些优化。

优化

首先可以把递归改为非递归版的迭代。其次一个重要优化是 蝴蝶变换 。他说的是我们递归的过程会把 $ F(x) $ 按奇偶分裂,我们最终分裂好的整个序列原来在 $ i $ (下标从 0 开始)位置的数,现在在 $ i $ 按照二进制翻转后的数为下标的地方,比如说 $ 1010111_2 (87_{10}) $ ,在翻转后是 $ 1110101_2 (117_{10}) $ 。

首先证明这个事:

我们脑动模拟一下递归过程,每次按照偶在左半边,奇在右半边去分,相当于看二进制下最后一位,按照大小确定了最终二进制下第一位的大小。

有点抽象,但是有图:

以下将经过 $ i $ 次分裂的序列称为 $ F_i $

这是一开始的 $ F_0(x) $

image

然后我们按照最后一位的大小,分开

image

这时可以看到 $ F_1 $ 的下标的最后一位和 $ F_0 $ 下标的第一位相同,因为一开始我们是按第一位大小排的,后来我们按最后一位大小排的。同时前一半和后一半没有关系了,他们各自内部都是按大小排好的,而且此时所有数的最后一位是啥不重要了,相当于所有数最后一位没了,倒数第二位顺延为最后一位,继续递归处理。

那么下一层也就根据 $ F_1 $ 的最后一位,也就是倒数第二位确定了第二位的是啥,然后 $ \cdots $ ,也就相当于二进制翻转。

image

而事实也确实如此。

好,设 $ pos_i $ 表示 $ i $ 翻转后的答案,那么我们怎么求 $ pos_i $ ? 可以考虑先翻转最后一位之前的,由 $ pos_{i/2} $ 得到,同时翻转最后最后一位,于是有 pos[i]=(pos[i>>1]>>1)|(i&1?n>>1:0) , $ n $ 是序列长度。

所以我们就获得了一个常数较小的 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 推荐更新“三次变两次优化”

三次变两次

因为是复数运算,设一个复数多项式 $ F(x) = f(x) + g(x) i $ ,那么有

\[F^2(x) = f^2(x) - g^2(x) + 2 f(x) g(x) i \]

所以我们要求的就是 $ F^2(x) $ 的虚部除以2,那么再插值的时候只需要插一次就行,从三次变成了两次,快了约 $ \frac{1}{4} $ 。

码(最终版)
#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] 力

给出序列 $ q $ ,有

\[E_i = \sum_{i=1}^{j-1} \frac{q_i}{(i-j)^2} - \sum_{i=j+1}^{n} \frac{q_i}{(i-j)^2} \]

求所有的 $ E_i $

我们考虑如何把 $ E_i $ 转成卷积的形式。首先把烦人的除法去掉,令 $ G(x) = \frac{1}{x^2} $ ,特殊的, $ G(0) = 0 $ ,同时令 $ q_0 = 0 $ ,那么有

\[E_i = \sum_{i=0}^{j} q(i) G(j-i) - \sum_{i=j}^{n} q(i) G(i-j) \]

前面已经是卷积的形式了,无需转化,后面的不是,我们考虑给他转成卷积的形式。

画个图看看:

image

他们是这样相乘的,但是卷积是这样相乘的:

image

所以我们考虑给 $ q(x) $ 转过来,变成 $ Q(x) $ ,所以后面的式子就变成了

\[\sum_{i=0}^{k=n-j} Q(i) G(k-i) \]

对两者 FFT 然后加起来即可。

二、 NTT

NTT ,中文 快速数论变换 是用来解决卷积过程中需要取模的快速加法卷积,说白了就是可以取模的 FFT。

这个东西就是直接把 $ w_n^1 $ 替换成了 $ g^{ \frac{mod-1}{n} } $ ,其中 $ mod $ 是模数,$ g $ 是模数的一个原根,$ n $ 是序列长度 ,但是原根这个东西我也不会,目前只知道 998244353 的原根是 3 ,所以 NTT 到此为止。

#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是用来解决有关位运算卷积的快速卷积算法。

即对于 $ c_i = \sum_{i = j \oplus k } a_j b_k $ ,在知道 $ a,b $ 的情况下,快速计算 $ c $ 。

或运算

此时 $ \oplus \to | $ 。

我们设 $ A,B,C $ 分别对应 $ a,b,c $ 进行了 FWT 后的序列,我们定义 $ A_i = \sum_{i=i|j} a_j $ ,那么自然有 $ C_i = A_i \cdot B_i $ 。

考虑如何快速求 $ A $ ,我们考虑初始序列 $ a $ ,它的长度 $ n $ 是 2 的整数次幂,我们将 $ a $ 一分为二,左边为 $ a_0 $ ,右边为 $ a_1 $ ,两者没有任何关系 ,直接递归求出 $ A_0 , A_1 $ ,此时需要需要合并 $ A_0,A_1 $ 为 $ A $ ,即,把 $ a_0 $ 向 $ A_1 $ 的贡献考虑进来, $ a_1 $ 向 $ A_0 $ 中的贡献考虑进来。

注意到将 $ a_1 $ 向 $ A_0 $ 考虑时,其中 $ a_1 $ 的目前最高位上一定是 1 ,而 $ a_0 $ 中的数一定是 0 ,所以 $ a_1 $ 不可能向 $ A_0 $ 做贡献。

再看 $ a_0 $ 向 $ A_1 $ 做贡献,注意到 $ \large a_{0_j} | a_{1_i} = a_{1_i} $ 当且仅当 $ \large a_{0_j} | a_{0_i} = a_{0_i} $ ,即能向 $ a_{0_i} $ 做贡献的都能向 $ a_{1_i} $ 做贡献,所以有

\[A = he( A_0 , A_0 + A_1 ) \]

其中 $ + $ 表示多项式相加, $ he(x,y) $ 表示将 $ x,y $ 两个序列顺序连接。所以我们可以分治去求 $ A $ 。

类似的可以退出逆运算为

\[A = he( A_0 , A_0 - A_1 ) \]

也就可以在得到 $ C $ 之后反推出 $ c $ 。


感觉就和 FFT , NTT 一样, FWT(|) 找了一种特殊运算,使得可以分治 $ O(n \log n) $ 的去求值插值。同时分治的过程可以看成逐渐考虑二进制上每一位的影响。


与运算

类比或运算,这里直接给出式子

求值:

\[A = he( A_0 + A_1 , A_1 ) \]

插值:

\[A = he( A_0 - A_1 , A_1 ) \]


异或运算

异或运算是 FWT 中最难理解的一个,是一种非常神奇的构造。

以下 $ \oplus \to xor $

首先我们定义 $ x \circ y = popcount( x \oplus y ) \mod 2 $ ,然后我们先给出 $ ( x \circ y ) \oplus ( x \circ z ) = x \circ ( y \oplus z ) $ 。

证明:

我们一位一位的去考虑,在一开始,两边都是 0,然后考虑新的一位,如果要做出改变,首先 $ x $ 在这位上得是 1,否则不用考虑。然后看 $ y \oplus z $ 的值,如果为 1 ,说明两个数一个是 1,一个是 0,所以此时左边右边的值都会改变,不影响等式成立。如果 $ y \oplus z = 0 $ ,说明两者一样,此时两边的值都不会改变,也不影响等式成立。

综上, $ ( x \circ y ) \oplus ( x \circ z ) = x \circ ( y \oplus z ) $ 。

接下来我们定义 $ A_i = \sum_{i \circ j = 0} a_j - \sum_{i \circ j = 1} a_j $ ,那么有

\[\begin{aligned} A_i \cdot B_i &= ( \sum_{i \circ j = 0} a_j - \sum_{i \circ j = 1} a_j ) \cdot ( \sum_{i \circ j = 0} b_j - \sum_{i \circ j = 1} b_j ) \\ &= (\sum_{i \circ j = 0} a_j \sum_{i \circ k =0} b_k + \sum_{i \circ j = 1} a_j \sum_{i \circ k =1} b_k) - ( \sum_{i \circ j = 0} a_j \sum_{i \circ k =1} b_k + \sum_{i \circ j = 1} a_j \sum_{i \circ k =0} b_k ) \\ &= \sum_{ i \circ (j \oplus k) =0 } a_j b_k - \sum_{ i \circ (j \oplus k) =1 } a_j b_k \\ &= C_i \end{aligned} \]

接下来考虑如何快速计算 $ A $ ,我们仍然尝试将它分为 $ A_0,A_1 $ 递归求解,接下来考虑合并。

首先举一个例子:

假设 $ A $ 长度为 8,我们写成二进制形式

image

然后我们一分为二,递归求解,但是注意递归到下一层后,他们在本层的第一位将不再考虑,也就是

image

这也就是为什么之前里面说分治的过程是逐渐考虑二进制位的影响。

理解了这里,下面就很好理解了。

对于 $ A $ ,分为 $ A_0 ,A_1 $ 之后,我们新考虑一位,定义 $ x' $ 为 考虑新一位的 $ x $ ,那么有 $ a_{0_i} ' \circ a_{1_j} ' = a_{1_i} \circ a_{1_j} , a_{0_i} ' \circ a_{0_j} ' = a_{0_i} \circ a_{0_j} , a_{1_i} ' \circ a_{1_j} ' = (a_{1_i} \circ a_{1_j}) \oplus 1 , a_{1_i} ' \circ a_{0_j} ' = a_{0_i} \circ a_{0_j} $ ,于是我们有

\[A = he( A_0 + A_1 , A_0 - A_1 ) \]

同时有逆推

\[A = he( \frac{ A_0 + A_1 }{2} , \frac{ A_0 - A_1 }{2} ) \]

模板题

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<<' ';
}

小结

感觉三者都是通过构造找了一些特殊的点或特殊的运算,使得我们可以通过这写特殊的性质达到快速计算的效果,当然这一切的前提都是因为我们用到的 $ O(n) $ 点乘对点是没有要求的,这就支撑着我们可以通过构造特殊点来完成快速卷积。

UPD:感谢各位指出本篇的 Huge 量错误。

posted @ 2025-02-09 20:37  lzrG23  阅读(443)  评论(19编辑  收藏  举报