多项式相关

合集 - 学习笔记(12)

1.并查集2023-11-24

2.多项式相关2023-11-24

3.CDQ分治2023-12-174.数论2023-12-155.数据结构2024-01-276.树分治2024-09-237.斜率优化（李超树）2024-06-268.贪心2024-08-069.单调队列优化dp2024-10-3110.扫描线2024-10-2911.网络流&费用流&二分图2024-10-2912.深搜DFS2023-03-25

收起

鸽了，$L_aT^eX$ 暂时不想修

虽然数学或许考得不多，但是考了不会就完了。

多项式快速乘法

当我们想对两个多项式做卷积的时候，朴素算法就是一个一个相乘，时间复杂度 $O(n^2)$。
但是有可以做到 $O(n\log n)$ 完成求卷积操作的算法，这里记 FFT 和 NTT ，这两种算法本质上只是快速地求出 DFT 和 IDFT 而已。
很显然我并不能完完全全的理解 DFT/IDFT 利用单位根运算的过程，所以应该背代码就行了，我不信会有人在这上面搞事情。

多项式的点值表示法

给定一个不超过 $n$ 次的多项式 $A(n)$ 以及 $n+1$ 个不同的点 $x_0,\ldots,x_n$，令 $y_i=A(x_i)$，则这 $n+1$ 组 $(x_i,y_i)$ 唯一的确定了这个多项式 $A(n)$。
如果我们得到两个多项式 $A(x),B(x)$ 分别在 $x_0,\ldots,x_n$ 下的点值，那么我们就可以在 $O(n)$ 的时间内求出 $A(x)\pm B(x)$ 或 $A(x)B(x)$ 在同一组位置的点值。
现在我们的问题是如何在多项式的系数和点值表示之间转换，这促使我们考虑一组特殊的点值。

补充内容：复数

自己去看高中数学，反正这里我们要用到他来计算。

单位根

这几张图片只是一些基本定义和性质，暂时和FFT没关系。

离散傅里叶变换（DFT）

把系数表示转为点值表示的过程。
定义：设 $a$ 是长度为 $n$ 的数列，对 $0\le k<n$，令 $b_k=\sum\limits_{i=0}^{n-1} a_i·\omega_n^{ki}$，则称 $b$ 为 $a$ 的离散傅里叶变换（DFT）。
可以看出，令 $A(x)=\sum a_ix^i$，则 $b_k$ 就是 $A(x)$ 在 $\omega_n^k$ 处的点值。因此计算 $a$ 的DFT与计算 $A(x)$ 在 $\omega_n^0,\omega_n^1,\ldots,\omega_n^{n-1}$ 处的点值表示是等价的。

离散傅里叶逆变换（IDFT）

这部分直接抄过来算了，知道是干什么的就行了。

目前为止，我们已经有了能将卷积这一部分的时间复杂度从 $O(n^2)$ 降至 $O(n)$ 的算法，但是DFT和IDFT的复杂度却是 $O(n^2)$ 的，这个时候我们就通过单位根来完成这两步操作。

FFT

FFT 算法的基本思想是分治。就 DFT 来说，它分治地来求当 $x=\omega_n^k$ 的时候 $f(x)$ 的值。基于FFT的分治思想体现在将多项式分为奇次项和偶次项处理。
idea is shit,show me your code：

void FFT(ll k,complex <db> *a,ll f){ //f表示我们做的是DFT还是IDFT
	if(k==1) return;
	complex <db> x[k/2+5],y[k/2+5];
	for(int i=0;i<=k;i+=2)
		x[i/2]=a[i],y[i/2]=a[i+1]; //分奇偶次项
	FFT(k/2,x,f);
	FFT(k/2,y,f);
	complex <db> W(cos(2.0*pi/k),f*sin(2.0*pi/k)),w(1.0,0.0); //切忌加等号
	for(int i=0;i<k/2;++i,w=w*W){ //使用单位根运算，背
		a[i]=x[i]+w*y[i];
		a[i+k/2]=x[i]-w*y[i];
	}
}

递推实现暂时不想写。
考虑到分治 DFT 能处理的多项式长度只能是 $2^m(m \in \mathbf{N}^ \ast )$，否则在分治的时候左右不一样长，右边就取不到系数了。所以要在第一次DFT之前就把序列向上补成长度为 $2^m(m \in \mathbf{N}^\ast )$（高次系数补 $0$）、最高项次数为 $2^{m-1}$ 的多项式。

NTT

大致流程基本和FFT一样，不过我们不再使用单位根来完成快速求DFT/IDFT，我们将使用一种神秘的东西：原根（不是原神）。

大概也是直接背代码就行了：

void NTT(ll k,ll *a,ll f){
	if(k==1) return;
	ll x[k/2+5],y[k/2+5];
	for(int i=0;i<=k;i+=2)
		x[i/2]=a[i],y[i/2]=a[i+1];
	NTT(k>>1,x,f),NTT(k>>1,y,f);
	ll w=1,W=(f==1?qpow(g,(mod-1)/k):qpow(inv,(mod-1)/k));
	for(int i=0;i<k/2;++i,w=w*W%mod){
		a[i]=(x[i]+w*y[i]%mod)%mod;
		a[i+k/2]=(x[i]-w*y[i]%mod+mod)%mod;
	}
}

还有很玄学的是 $k$ 的取值，以后尽量这么来：ll base=ceil(log2(n+m+1)); ll k=1<<base;。直接while来不断乘二会莫名RE，但是FFT就不会？

总结

这些东西是绝对不可能单独考的，肯定要结合其他数学知识，做题的时候重在将题目中的式子转化成我们可以进行快速计算的形式。
还有些诡异的东西比如任意模数多项式乘法暂时不想搞了，数学内容慢慢来，虽然似乎没有多少时间让我慢慢来了……

如何推式子？

不知道
我们一般要尽量把题目中给定的式子转化成卷积形式：$C_k=\sum\limits_{i=0}^k A_i B_{k-i}$，题目中的式子一般都长得非常神秘，我们可以考虑从以下方面入手。

• 式子中有多个求和考虑把它们分别都化成卷积形式

• 改变式子中 $sum$ 的上下标，比如把不等式改成等价的等式，还有注意求和区间的开闭情况，最好是写成闭区间的形式。

• 找到式子中的不变项把它提出来，或是把能消去的部分都给消去，把式子化得越简单越好。

• 设元代替式子中的某些单项式部分，然后观察这样的话更容易发现怎么转化。

• 某些部分珂以考虑把它展开后分析

• 当我们化出来后某一部分的下标变成了形如 $i+j$ 的形式，珂以考虑反转序列额外开一个数组记录，比如我们原来是 $\sum\limits_{i=0}^n f[i+j]*g[i]$，珂以令 $t=n-j$，并且 $fx$ 为 $f$ 的反转序列，然后转化成 $\sum\limits_{i=0}^n fx[t-i]*g[i]$

• 珂以将一个抽象的式子转化成别的东西，比如组合数什么的，然后再化成卷积形式

• 有些题可以考虑将每项的贡献计算出来，或是被计算次数，然后与原数列直接卷积。

例题

P3338 [ZJOI2014] 力

滥用 $LaTeX$ 警告
首先我们直接写成 $E_j=\sum\limits_{i=\small{1}}^{j-\small{1}}\dfrac{q_j}{(i-j)^{\small{2}}}-\sum\limits_{i=j+\small{1}}^{n}\dfrac{q_j}{(i-j)^{\small{2}}}$，其中我们把 $q_i$ 项消去了。
然后我们令 $f_i=q_i,g_i=\dfrac{\small{1}}{i^{\small{2}}}$，则式子珂以写成 $\sum\limits_{i=\small{1}}^{j-\small{1}} f[i]*g[j-i]-\sum\limits_{i=j+\small{1}}^n f[i]*g[j-i]$
考虑改变 $\sum$ 的上下标，令 $f[\small{0}]=g[\small{0}]=\small{0}$，同时我们能把开区间换成闭区间（虽然按式子的形式来说这样子就会除以 $\small{0}$，但是对于实际计算是不影响的）那么又珂以写成 $\sum\limits_{i=\small{0}}^{j} f[i]*g[j-i]-\sum\limits_{i=j}^n f[i]*g[j-i]$
好的，我们已经将左边的求和化成了卷积形式，现在只关注右边。
我们将右边的式子展开：$f[i]*g[\small{0}]+f[i+\small{2}]*g[\small{1}]+\ldots+f[i+(n-j)]*g[n-j]$，
那么考虑将式子写成 $\sum\limits_{i=\small{0}}^{n-j} f[j+i]*g[i]$
好，然后我们用到上面所述的反转序列的技巧，令 $fx$ 为 $f$ 的反转序列，同时 $t=n-j$，
这样就可以转化成 $\sum\limits_{i=\small{0}}^t fx[t-i]*g[i]$
好！这样子最终的式子就是这样的：$\sum\limits_{i=\small{0}}^{j} f[i]*g[j-i]-\sum\limits_{i=\small{0}}^t fx[t-i]*g[i]$，已经是一个可以拿来FFT的东西了。
我们再设 $a[i]=f[i],b[i]=fx[i],c[i]=g[i]$ 这样把这三个序列经过DFT后的运算就是 $a_i=a_i*c_i,b_i=b_i*c_i$，最后输出 $b_i-a_i$ 就行了。

code

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define db double
const ll N=3*114514,M=1919810;
const db pi=3.1415926535897932384626;
ll n; db w;
complex <db> a[N],b[N],c[N];
void FFT(ll k,complex <db> *a,ll f){
	if(k==1) return;
	complex <db> x[k/2+5],y[k/2+5];
	for(int i=0;i<=k;i+=2)
		x[i/2]=a[i],y[i/2]=a[i+1];
	FFT(k>>1,x,f),FFT(k>>1,y,f);
	complex <db> W(cos(2.0*pi/k),f*sin(2.0*pi/k)),w(1.0,0.0);
	for(int i=0;i<k/2;++i,w=w*W){
		a[i]=x[i]+w*y[i];
		a[i+k/2]=x[i]-w*y[i];
	}
}
int main(){
	cin>>n;
	for(int i=1;i<=n;++i){
		cin>>w;
		a[n-i].real(w),b[i].real(w);
		c[i].real((1.0/i*1.0)*(1.0/i*1.0));
	}
	ll k=1;
	while(k<=2*n) k<<=1;
	FFT(k,a,1),FFT(k,b,1),FFT(k,c,1);
	for(int i=0;i<=k;++i) a[i]=a[i]*c[i],b[i]=b[i]*c[i]; //点值表示是0~k 
	FFT(k,a,-1),FFT(k,b,-1);
	for(int i=1;i<=n;++i) printf("%0.10lf\n",b[i].real()/k-a[n-i].real()/k);
    //记得递归版的FFT要对每一项除以k
	return 0;
}

终于写完了，但是还有

P5641 【CSGRound2】开拓者的卓识

卡老师出的好题，难。
考虑贡献法计算每个 $a_i$ 对 $sum$ 的贡献，也就是我们要找到有多少个一阶子段和包含 $a_i$。然后组合数就推不出来了……
考虑怎么去找这些一阶子段和，也就是区间个数。我们可以理解为是去找左右端点，满足 $1\le l_{\small{1}}\le l_{\small{2}}\le\ldots\le l_k\le i,i\le r_k\le r_{k-\small{1}}\le\ldots\le r_{\small{1}}$。
那么对于总的答案就是要去找满足 $1\le l_{\small{1}}\le l_{\small{2}}\le\ldots\le l_k\le n$ 的 $(l_{\small{1}},l_{\small{2}},\ldots,l_k)$。
妈的推不来组合数，反正最后的结论是 $sum_{k,\small{1},r}=\sum\limits_{i=\small{1}}^n a_i\times C\binom{i+k-\small{2}}{k-\small{1}}C\binom{r-i+k-\small{1}}{k-\small{1}}$，然后考虑怎么把它化成卷积形式。
这个就比较简单了，我们令 $A_i=a_{i+\small{1}}\times C\binom{i+k-\small{1}}{k-\small{1}},B_i=C\binom{i+k-\small{1}}{k-\small{1}},i\in[\small{0},r)$，然后答案就是 $\sum\limits_{i=\small{0}}^{r-\small{1}}A_i*b_{r-i-\small{1}}$，已经是个卷积的形式了。
然后就预处理一下 $C\binom{i+k-\small{1}}{k-\small{1}}$ 的值就行了，有取模所以用NTT，注意有取模的时候不能有÷而是乘逆元（这常识我竟然现在才知道/qd）。

code

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll N=114514,M=4*1919810;
const ll mod=998244353,g=3,inv=332748118;
ll n,m,a[M],b[M];
ll qpow(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
} 
void NTT(ll k,ll *a,ll f){
	if(k==1) return;
	ll x[k/2+5],y[k/2+5];
	for(int i=0;i<=k;i+=2)
		x[i/2]=a[i],y[i/2]=a[i+1];
	NTT(k>>1,x,f),NTT(k>>1,y,f);
	ll w=1,W=(f==1?qpow(g,(mod-1)/k):qpow(inv,(mod-1)/k));
	for(int i=0;i<k/2;++i,w=w*W%mod){
		a[i]=(x[i]+w*y[i]%mod)%mod;
		a[i+k/2]=(x[i]-w*y[i]%mod+mod)%mod;
	}
}
int main(){
	cin>>n>>m;
	for(int i=0;i<n;++i) cin>>a[i];
	ll ba=ceil(log2(2*n+1)),k=1<<ba;
	b[0]=1;
	for(int i=1;i<n;++i) b[i]=(b[i-1]*(i+m-1)%mod)*qpow(i,mod-2)%mod;//逝了,果然还是要乘逆元 
	for(int i=0;i<n;++i) a[i]=a[i]*b[i]%mod;
	NTT(k,a,1),NTT(k,b,1);
	for(int i=0;i<k;++i) a[i]=a[i]*b[i]%mod;
	NTT(k,a,-1);
	ll invv=qpow(k,mod-2);
	for(int i=0;i<n;++i) cout<<a[i]*invv%mod<<" ";
	return 0;
}

P5488 差分与前缀和

其实不难，没上一道题难，还有这个 $k$ 这么大确实很鬼畜。
$k$ 阶差分/前缀和本质上就是求 $k$ 遍差分或前缀和而已，那么我们自然就想到处理出每一位会被算到的次数。
我不会严谨的证明，只会感性理解：对于一个数列 $a$，我们对它求 $k$ 次前缀和，那么珂以设 $b_i$ 为到 $i$ 位置时每个数被算了多少次，相当于我们同时也在对计算次数做前缀和。那么我们就珂以直接将得出的 $b$ 与 $a$ 做卷积，就可以得出每一项最终的答案。对于差分，因为差分是前缀和的逆运算，所以我们也珂以将 $b$ 数组倒过来求。

code

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll N=3*114514,M=1919810;
const ll mod=1004535809,g=3;
ll n,t,a[N],b[N];
string k_; //鬼畜
ll qpow(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}
void NTT(ll k,ll *a,ll f){
	if(k==1) return;
	ll x[k/2+1],y[k/2+1];
	for(int i=0;i<=k;i+=2)
		x[i/2]=a[i],y[i/2]=a[i+1];
	NTT(k>>1,x,f),NTT(k>>1,y,f);
	ll w=1,W=qpow(g,(mod-1)/k);
	if(f==-1) W=qpow(W,mod-2)%mod;
	for(int i=0;i<k/2;++i,w=w*W%mod){
		a[i]=(x[i]+w*y[i]%mod)%mod; //写错了/fn 
		a[i+k/2]=(x[i]-w*y[i]%mod+mod)%mod;
	}
}
int main(){
	ios::sync_with_stdio(0);
	cin.tie(0); cout.tie(0);
	cin>>n>>k_>>t; --n;
	for(int i=0;i<=n;++i) cin>>a[i];
	ll ba=ceil(log2(2*n+1)),k=1<<ba;
	ll id=1,K=0;
	for(int i=k_.size()-1;i>=0;--i)
		K=(K+(k_[i]-'0')*id)%mod,id=id*10%mod;
	b[0]=1;
	if(!t) for(int i=1;i<=n;++i) b[i]=b[i-1]*qpow(i,mod-2)%mod*(K+i-1)%mod;
	else for(int i=1;i<=n;++i) b[i]=(-b[i-1]*qpow(i,mod-2)%mod*(K-(i-1))%mod+mod)%mod;
	//for(int i=0;i<=n;++i) cout<<b[i]<<" ";cout<<'\n';
	NTT(k,a,1),NTT(k,b,1);
	for(int i=0;i<=k;++i) a[i]=a[i]*b[i]%mod;
	NTT(k,a,-1);
	ll inv=qpow(k,mod-2);
	for(int i=0;i<=n;++i) cout<<a[i]*inv%mod<<" ";
	return 0;
}

生成函数

生成函数是一种用来描述特定组合对象的形式多项式，组合对象满足特定性质，如大小长度为 $n$ 等，其描述的组合对象数是有限的。其组合对象通常分为有标号和无标号两类，无标号主要描述的就是组合，有标号就描述排列。关于一般生成函数的组合式意义：对于一个数列在其中选数（不考虑顺序）第 $n$ 项的系数表示方案数，每项次数表示选取的个数。如果两数列相加便是生成函数（多项式）的相加，两数列的卷积便是生成函数的乘积。

普通生成函数（OGF）

描述数列 $A_0,A_1,\ldots$ 的普通生成函数（OGF）定义为形式幂级数 $A(x)=\sum\limits_{i=0}^{\infty} A_i x^i$，将其写成级数形式为 $a_0x^0+a_1x^1+a_2x^2+\ldots$，根据等比数列公式我们珂以将其写成封闭形式为 $\frac{a_i}{(a_i-x)^i}$。这种生成函数用来描述无标号对象。

指数型生成函数（EGF）

描述数列 $A_0,A_1,\ldots$ 的指数型生成函数（EGF）定义为形式幂级数 $A(x)=\sum\limits_{i=0}^{\infty}\frac{A_ix^i}{i!}$，用于描述有标号对象。考虑有标号组合对象的拼接，将对象 $a,b$ 拼接为 $c$，其中 $size(a)=n,size(b)=m$，那么由于 $a,b$ 间的相对顺序会变化，于是我们会有 $\binom{n}{m}=\frac{(n+m)!}{n!m!}$ 种排列方法，因此若将组合对象 $A,B$ 拼接为 $C$，则有 $C_n=\sum\limits_{i+j=n}A_iB_j\frac{(i+j)!}{i!j!}$，则有 $C(x)=A(x)B(x)$

例题

P2000 拯救世界

根据题意，得到一下十个生成函数：

把他们全部乘起来，得到了一个很简单的生成函数 $\frac{1}{(1-x)^5}$，我们知道 $\frac{1}{(1-x)^n}=C^{i}_{n+i-1}x^i$，则这个式子就可以化成 $\frac{1}{(1-x)^5}=\sum\limits_{i=0}^\infty C_{i+4}^i x^i$，而系数 $C_{i+4}^i=\frac{(n+4)(n+3)(n+2)(n+1)}{4!}$，所以答案其实就是 $(n+\small{1})(n+\small{2})(n+\small{3})(n+\small{4})/ \small{24}$ 但是这个题的 $n$ 非常鬼畜，所以还要用FFT/NTT解决，注意我们不能直接一次点乘乘完而是要分四次，中间处理进位。这样子就总共进行了11次NTT，常数巨大，只能用递推版的，递归般的TLE了/oh

code

点击查看代码

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll N=8*114514,M=1919810;
const ll mod=998244353,g=3,G=332748118;
ll qpow(ll a,ll b){
	ll ans=1;
	while(b){
		if(b&1) ans=ans*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return ans;
}
ll r[N];
void NTT(ll n,ll *a,ll f){ //论常数 
	for(int i=0;i<n;++i)
		if(i<r[i]) swap(a[i],a[r[i]]);
	for(int i=1;i<n;i<<=1){
		ll w=qpow( f==1 ? g : G,(mod-1)/(i<<1));
		for(int j=0;j<n;j+=(i<<1)){
			ll W=1;
			for(int k=0;k<i;++k,W=W*w%mod){
				ll x=a[j+k],y=W*a[i+j+k]%mod;
				a[j+k]=(x+y)%mod;
				a[i+j+k]=(x-y+mod)%mod;
			}
		}
	}
}
string s;
ll n,a[N],b[N],c[N],d[N],ans[N],res[N];
int main(){
	cin>>s; n=s.size()-1;
	for(int i=1;i<=n;++i) a[i]=b[i]=c[i]=d[i]=s[n-i]-'0';
	a[0]=s[n]-'0'+1,b[0]=s[n]-'0'+2,c[0]=s[n]-'0'+3,d[0]=s[n]-'0'+4;++n;
	ll ba=ceil(log2(4*n+1)),k=1<<ba;
	ll inv=qpow(k,mod-2);
	for(int i=0;i<k;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(ba-1));
	NTT(k,a,1),NTT(k,b,1),NTT(k,c,1),NTT(k,d,1);
	for(int i=0;i<=k;++i) ans[i]=a[i]%mod;
	NTT(k,ans,-1);
	for(int i=0;i<=k;++i) ans[i]=ans[i]*inv%mod;
	for(int i=0;i<=k;++i) ans[i+1]+=ans[i]/10,ans[i]%=10;
	NTT(k,ans,1);
	for(int i=0;i<=k;++i) ans[i]=ans[i]*b[i]%mod;
	NTT(k,ans,-1);
	for(int i=0;i<=k;++i) ans[i]=ans[i]*inv%mod;
	for(int i=0;i<=k;++i) ans[i+1]+=ans[i]/10,ans[i]%=10;
	NTT(k,ans,1);
	for(int i=0;i<=k;++i) ans[i]=ans[i]*c[i]%mod;
	NTT(k,ans,-1);
	for(int i=0;i<=k;++i) ans[i]=ans[i]*inv%mod;
	for(int i=0;i<=k;++i) ans[i+1]+=ans[i]/10,ans[i]%=10;
	NTT(k,ans,1);
	for(int i=0;i<=k;++i) ans[i]=ans[i]*d[i]%mod;
	NTT(k,ans,-1);
	n*=4;
	for(int i=0;i<=n;++i) ans[i]=ans[i]*inv%mod;
	for(int i=0;i<n;++i) ans[i+1]+=ans[i]/10,ans[i]%=10;
	ll len=0;
	for(int i=n;i>=0;--i)
		if(ans[i]!=0){
			len=i;
			break;
		}
	for(int i=len;i>=0;--i){
		res[i]=ans[i]/24;
		ans[i]%=24;
		ans[i-1]+=ans[i]*10;
	}
	while(len&&res[len]==0) --len;
	for(int i=len;i>=0;--i) cout<<res[i];
	return 0;
} 

快速沃尔什变换（FWT）

OI中FWT是用于解决对下表进行位运算卷积问题的方法。$c_i=\sum\limits_{i=j\oplus k}a_jb_k$，$\oplus$ 表示二元位运算的一种，这里记录或、与、异或三种运算。

核心思想和FFT/NTT基本一致：我们令序列 $a$ 进行快速沃尔什变换后的序列为 $fwt[a]$，已知序列 $a,b$ 求新序列 $c=a·b$，直接计算 $O(n^{2})$，但我们珂以先 $O(n\log n)$ 求出 $fwt[a],fwt[b]$ 再 $O(n)$ 求出 $fwt[c]=fwt[a]·fwt[b]$ 再 $O(n\log n)$ 逆变换求出 $c$，总时间复杂度为 $O(n\log n)$。

或&与

要求 $c_i=\sum\limits_{i=j|k} a_jb_k$，可以得到若 $j|i=i,k|i=i$ 则 $(j|k)|i=i$，考虑构造 $fwt[a]_i=\sum_{j|i=i}a_j$，则有 $fwt[a]\times fwt[b]=\sum\limits_{j|i=i}\sum\limits_{k|i=i}a_jb_k$

​
$\begin{alignat*}{3} fwt[a]\times fwt[b]&=\sum\limits_{j|i=i}\sum\limits_{k|i=i}a_jb_k\\ &=\sum\limits_{(j|k)|i=i}a_jb_k\\ &=fwt[c] \end{alignat*}$
现在我们的问题是怎么求出 $fwt[a]$，我们令 $a_{\small{0}}$ 表示 $a$ 中下标最高位为 0 的那部分序列，$a_{\small{1}}$ 表示下标最高位为 1 的那部分序列。

则有：$fwt[a]=\mathrm{merge}(fwt[a_{\small{0}}],fwt[a_{\small{0}}]+fwt[a_{\small{1}}])$，其中 merge 为拼接两个数列，$+$ 表示对应位置相加。

于是我们珂以类似FFT/NTT的分治方式处理：

void Or(ll *a,ll n,ll f){
	for(int m=2,k=1;m<=n;m<<=1,k<<=1)
		for(int i=0;i<m;i+=m)
			for(int j=0;j<k;++j)
				a[i+j+k]=((a[i+j+k]+f*a[i+j])%mod+mod)%mod;
}

与的推导过程和运算与之同理：

void And(ll *a,ll n,ll f){
	for(int m=2,k=1;m<=n;m<<=1,k<<=1)
		for(int i=0;i<m;i+=m)
			for(int j=0;j<k;++j)
				a[i+j]=((a[i+j]+f*a[i+j+k])%mod+mod)%mod;
}

异或

我们定义 $x\otimes y=\text{popcount}(x\&y)\pmod 2$，$\text{popcount}$ 表示二进制下 $\mathtt 1$ 的个数。

满足 $(i\otimes j)xor(i\otimes k)=i\otimes(j\:xor\:k)$

构造 $fwt[a]_i=\sum\limits_{i\otimes j=\mathtt{0}}a_j-\sum\limits_{i\otimes j=\mathtt{1}}a_j$

则有：（太长了放图片）

因此有
$\begin{alignat*}{2} fwt[a]&=\text{merge}(fwt[a_0]+fwt[a_1],fwt[a_0]-fwt[a_1])\\a&=\text{merge}(\frac{a_0+a_1}{2},\frac{a_0-a_1}{2})\end{alignat*}$
注意逆变换的时候乘的是 $2$ 的逆元。
代码：

void Xor(ll *a,ll n,ll f){
	for(int m=2,k=1;m<=n;m<<=1,k<<=1)
		for(int i=0;i<m;i+=m)
			for(int j=0;j<k;++j){
				a[i+j]=(a[i+h]+a[i+j+k])%mod;
				a[i+j+k]=(a[i+j]-a[i+j+k]*2%mod+mod)%mod;
				a[i+j]=a[i+j]*f%mod,a[i+j+k]=a[i+j+k]*f%mod;
			}
}