FFT/FWT 相关理论推导复习

本文下标一般从 $0$ 开始，部分除法用 $/$。

建议前置知识：矩阵乘法，复数四则运算，单位根，模素数下原根。其实也完全可以在阅读过程中到不会的当场搜索即可，需要有点矩阵乘法基础知识。

本篇文章完成周期为一周，每个式子均为作者从头到尾手推过的。所以一下子看不懂可以别急，多看看。代码中均有注释，实在看不懂可以对着代码理解。

说实话可能这适合看过一两遍 FFT 的人再来学，但是可能之后我介绍 FFT 就用这篇了。

卷积：记长为 $n$ 的数组 $a,b$ 在运算 $\circ$ 下的卷积 $a\circ b=c$，其中 $c_k=\sum\limits_{i\circ j=k} a_ib_j$，$c$ 也为长为 $n$ 的数组。

直接暴力计算卷积复杂度为 $O(n^2)$，其中 $n$ 为数组长度。

DFT-IDFT#

一般快速计算特殊卷积的方法为构造 DFT 变换：

欲构造可逆的 DFT 变换 把长度为 $n$ 的数组映射到长度为 $n$ 的数组，满足 $\texttt{DFT}(a)_i\times \texttt{DFT}(b)_i=\texttt{DFT}(c)_i$。

即把卷积通过 DFT 变换 转化成对应位置的点乘。

设 DFT 的逆变换 $\texttt{IDFT}=\texttt{DFT}^{-1}$，此时 $c=\texttt{IDFT}(\{\texttt{DFT}(a)_i\times \texttt{DFT}(b)_i\})$，于是我们只需快速做 $\texttt{DFT},\texttt{IDFT}$ 变换即可。

一般代码流程为：把 $a,b$ 做 DFT 变换 存到原来的 $a,b$ 上，对应位点乘：$\forall i,a_i\gets a_i\times b_i$，把 $a$ 做 IDFT 变换 后得到的结果就是 $c$ 了。

一般地，不妨设 DFT 变换为线性变换，于是设 $\vec{u}=\begin{bmatrix} a_0\\ a_1\\ \vdots\\ a_{n-1} \end{bmatrix}$，$v,w$ 为 $b,c$ 同理写成向量。

则 DFT 变换 即为把 $\vec{u}$ 左乘一个 $n\times n$ 矩阵：$\texttt{DFT}(\vec{u})=A\vec{u}$，其中 $A$ 为 $n\times n$ 的可逆矩阵。

此时需要满足 $(A\vec{u})_i\times (A\vec{v})_i=(A\vec{w})_i$。

做点一般性推导：$(A\vec{u})_i\times (A\vec{v})_i=(A\vec{w})_i\Leftrightarrow \sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}A_{i,j}A_{i,k} a_jb_k=\sum\limits_{j=0}^{n-1}\sum\limits_{k=0}^{n-1}A_{i,j\circ k} a_jb_k$。

于是 $A$ 需要满足：$\forall 0\le i,j,k<n,A_{i,j}A_{i,k}=A_{i,j\circ k}$。称此时的 $A$ 为 DFT 矩阵。

FFT/NTT#

快速计算 $a\times b=c$，其中 $c_k=\sum\limits_{i+j=k} a_ib_j$。模板题。

注：若令 $f(x)=\sum\limits_{i=0}^n a_ix^i,g(x)=\sum\limits_{i=0}^m b_ix^i,h(x)=\sum\limits_{i=0}^{n+m} c_ix^i$，则我们即需求出 $h(x)=f(x)g(x)$，即求多项式乘法。

范德蒙德矩阵#

给定 $x_0,x_1,\cdots,x_{n-1}$，有 $n\times n$ 的范德蒙德矩阵为：$A:A_{i,j}=x_{j}^i$，即 $A=\begin{bmatrix} 1&1&\cdots & 1\\ x_0&x_1&\cdots & x_{n-1}\\ \vdots&\vdots&\vdots&\vdots\\ x_0^{n-1}&x_1^{n-1}&\cdots & x_{n-1}^{n-1}\\ \end{bmatrix}$。

设 $w_n$ 为 $n$ 次单位根，则有性质：$w_{n}^n=1,w_{n}^{n/2}=-1,w_n^k=\cos\left(\dfrac{2k\pi}{n}\right)+\mathrm {i}\sin\left(\dfrac{2k\pi}{n}\right)$。

此时有人类智慧：令上文中 $x_i=w_n^i$，则此时范德蒙德矩阵 $A:A_{i,j}=w_n^{ij}$ 即为 DFT 变换 所需的矩阵。

证明：考虑一般性推导，此时 $A_{i,j}A_{i,k}=w_n^{ij}w_n^{ik}=w_n^{i(j+k)}=A_{i,j+k}$。

注意到此时 $c$ 的项数不对，设 $a,b$ 的最高次分别为 $n,m$，在 $i>n$ 时的 $c_i$ 依然需要计算。不妨暂时把 $a,b,c$ 都看成长度为 $n+m$ 的数组，不足补 $0$ 即可。

---

为何不可取 $A:A_{i,j}=w_n^{j}$，此时由上述证明显然也满足 $(A\vec{u})_i\times (A\vec{v})_i=(A\vec{w})_i$。因为此时 $A$ 的行两两相等，$\det A=0$，即 $A$ 不可逆，不满足条件。

考虑如何求逆矩阵：观察 $(A\times A)_{i,j}=\sum\limits_{k=0}^{n-1} w_{n}^{ik}w_{n}^{kj}=\sum\limits_{k=0}^{n-1} w_{n}^{(i+j)k}=\dfrac{1-w_n^{(i+j)n}}{1-w_n^{(i+j)}}$，注意到一般情况下分母为 $0$，分子非 $0$。$(\star)$

做改动：设矩阵 $B:B_{i,j}=w_n^{-ij}$，则类似的：$(A\times B)_{i,j}=\begin{cases}n(i=j)\\ \dfrac{1-w_n^{(i-j)n}}{1-w_n^{(i-j)}}=0(i\neq j)\end{cases}$

于是此时 $A^{-1}=\dfrac{1}{n}B$，$A\times A^{-1}$ 就被构造为了单位矩阵。

快速做 DFT/IDFT#

设 $w_n$ 为 $n$ 次单位根，$A:A_{i,j}=w_n^{ij}$，欲快速求 $\vec{U}=A\vec{u}$，此时 $U_i=\sum\limits_{j=0}^{n-1} w_n^{ij}a_j$。

人类智慧：把 $n\gets 2^k$ 补齐为 $2$ 的幂次，对这个过程分治：

把 $U_i=\sum\limits_{j=0}^{n-1} w_n^{ij}a_j$ 看做把 $x=w_n^{i}$ 带入多项式 $f(x)=\sum\limits_{j=0}^{n-1} a_jx^j$。

分治：设 $f(x)=f_0(x^2)+xf_1(x^2)$，其中 $f_0(x)=\sum\limits_{j=0}^{n/2-1} a_{2j}x^j,f_1(x)=\sum\limits_{j=0}^{n/2-1} a_{2j+1}x^j$。

递归下去，不妨设已经求得 $U_{0,i}=f_0(w_{n/2}^i),U_{1,i}=f_1(w_{n/2}^i)$。考虑求 $U$：

$U_i=f(w_n^{i})=f_0(w_n^{2i})+w_n^{i}f_1(w_n^{2i})=f_0(w_{n/2}^{i})+w_n^{i}f_1(w_{n/2}^{i})=U_{0,i}+w_n^{i}U_{1,i}$

$U_{i+n/2}=f(w_n^{i+n/2})=f_0(w_n^{2i+n})+w_n^{i+n/2}f_1(w_n^{2i+n})=f_0(w_{n/2}^{i})-w_n^{i}f_1(w_{n/2}^{i})=U_{0,i}-w_n^{i}U_{1,i}$

于是只需要在分治的过程中实现一下复数类（要求复数四则运算），然后按上述过程做即可。

注意下 $2^k$ 的选取，设看成多项式时，$a,b$ 的最高次数为 $N,M$，则之前取 $n=N+M$，于是应有 $2^k>N+M$。取不取等啥的细节取决于你的实现。

---

IDFT 同理做即可，$f,U$ 等记号类似同上，只不过 $U_i=\sum\limits_{j=0}^{n-1} w_n^{-ij}a_j$。

并且把 $U_i=\sum\limits_{j=0}^{n-1} w_n^{-ij}a_j$ 看做把 $x=w_n^{-i}$ 带入多项式 $f(x)=\sum\limits_{j=0}^{n-1} a_jx^j$。

$U_i=U_{0,i}+w_n^{-i}U_{1,i},U_{i+n/2}=U_{0,i}-w_n^{-i}U_{1,i}$。

由于 $A^{-1}=\dfrac{1}{n}B$，注意最后要 $/n$，而且除的是补齐后的 $2^k$ 并非原始的 $n$。复杂度 $O(n\log n)$。

但是暴力分治做常数太大了，于是有蝴蝶变换。稍后讲。

FFT 三次变两次优化#

欲求两个整系数多项式 $f,g$ 的卷积 $h=f\times g$。直接朴素做 DFT-IDFT 需要对 $f,g$ 分别做一次 DFT，对乘后的数组做一次 IDFT。总共要三次 DFT/IDFT。

考虑由于需要单位根，FFT 需要在复数类下实现。而一般需要做的卷积都是整系数多项式，于是可以充分利用复数性质。

考虑 $(f(x)+\mathrm {i}g(x))^2=f(x)^2-g(x)^2+2f(x)g(x)\mathrm {i}$，于是 $f(x)g(x)=\dfrac{1}{2}{\operatorname{Im}}((f(x)+\mathrm {i}g(x))^2)$。

于是把 $g$ 放在 $f$ 的虚部上，对 $f$ 做一次 DFT，对 $f$ 的每一项点值平方，对 $f$ 做一次 IDFT，取虚部每一项 $/2$ 即可。总共要两次 DFT/IDFT。

由于 FFT 常数较大，于是减少一次优化显著。同时由于 DFT/IDFT 的形式相似，于是代码中我写成了一个函数。

$\texttt{code}$

#include<bits/stdc++.h>
#define LL long long
#define LD double
using namespace std;
const int N=4e6+5;
const LD Pi=acos(-1.0);
int mmax,n,m,o;
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}//bger(n) 表示 >n 的第一个 2 的幂次
struct comp
{
	LD x,y;
	inline comp operator+(const comp &xx)const{return {x+xx.x,y+xx.y};}
	inline comp operator-(const comp &xx)const{return {x-xx.x,y-xx.y};}
	inline comp operator*(const comp &xx)const{return {x*xx.x-y*xx.y,x*xx.y+y*xx.x};}
}a[N],b[N],w[N];//写个复数类，支持加减乘 
void DFT(int l,int r)
{
	if(l==r) return;int mid=(l+r)>>1;
	for(int i=l;i<=r;i++) b[i]=a[i];int t=l-1;
	for(int i=l;i<=r;i+=2) a[++t]=b[i];
	for(int i=l+1;i<=r;i+=2) a[++t]=b[i];//提取奇偶项 
	DFT(l,mid);DFT(mid+1,r);int L=mid-l+1;//递归求解 
	for(int i=l,k=0;i<=mid;i++,k++)
	{
		comp u=a[i],w={cos(o*Pi*k/L),sin(o*Pi*k/L)},v=a[i+L]*w;
		a[i]=u+v,a[i+L]=u-v;//转移合并出现在的 U,即这里的 a,
		//注意这里 pi*k/L 实际上是 2pi*k/(2L),2L 才是当前长度 n 
	}
}
int main()
{
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i].x;
	for(int i=0;i<=m;i++) cin>>a[i].y;mmax=bger(n+m);
	o=1;DFT(0,mmax-1);
	for(int i=0;i<mmax;i++) a[i]=a[i]*a[i];
	o=-1;DFT(0,mmax-1);
	for(int i=0;i<=n+m;i++) cout<<(int)(a[i].y/2/mmax+0.5)<<" ";
	return 0;
}

蝴蝶变换#

可是上面的 FFT 还是有点慢，考虑不递归实现。

只考虑上面实现中 DFT 的前半部分，观察所有递归下去的奇偶项交换都执行完后的 $a$ 数组会如何变换。

void DFT(int l,int r)
{
	if(l==r) return;int mid=(l+r)>>1;
	for(int i=l;i<=r;i++) b[i]=a[i];int t=l-1;
	for(int i=l;i<=r;i+=2) a[++t]=b[i];
	for(int i=l+1;i<=r;i+=2) a[++t]=b[i];
	DFT(l,mid);DFT(mid+1,r);int L=mid-l+1;
}

记 $n$ 为补齐 $2$ 的幂次后的数组长度。观察 $n=8$ 的情况：

称这个「位逆序置换」为蝴蝶变换。记 $i$ 变换后的位置为 $rev_i$，则变换过程即为交换所有位置 $(i,rev_i)$。考虑如何求 $rev_{0,1,\cdots,n-1}$。

对着代码理解即可

// 同样需要保证 len 是 2 的幂
// 记 rev[i] 为 i 翻转后的值
void change(Complex y[], int len) {
  for (int i = 0; i < len; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) {  // 如果最后一位是 1，则翻转成 len/2
      rev[i] |= len >> 1;
    }
  }
  for (int i = 0; i < len; ++i) {
    if (i < rev[i]) {  // 保证每对数只翻转一次
      swap(y[i], y[rev[i]]);
    }
  }
  return;
}

之后需要枚举层数，然后合并出对应的 $i,i+n/2$ 位置。代码挺短的，看代码理解吧。

$\texttt{code}$

#include<bits/stdc++.h>
#define LL long long
#define LD double
using namespace std;
const int N=4e6+5;
const LD Pi=acos(-1.0);
int mmax,n,m,rev[N];
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}//bger(n) 表示 >n 的第一个 2 的幂次
struct comp
{
	LD x,y;
	inline comp operator+(const comp &xx)const{return {x+xx.x,y+xx.y};}
	inline comp operator-(const comp &xx)const{return {x-xx.x,y-xx.y};}
	inline comp operator*(const comp &xx)const{return {x*xx.x-y*xx.y,x*xx.y+y*xx.x};}
}a[N],b[N],w[N];//写个复数类，支持加减乘 
void DFT(comp *a,int n,int o)
{
	for(int i=1;i<n;i++) rev[i]=rev[i>>1]>>1,(i&1)&&(rev[i]|=n/2);//求 rev 
	for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);//交换，注意别交换两次了 
	for(int i=1;i<n;i<<=1)//枚举层 
	{
		for(int j=0;j<n;j+=(i<<1))//枚举区间起点 
		{
			for(int k=j,K=0;k<j+i;k++,K++)//枚举区间中每一个点
			//这里建议自己画一个 n=8 的图理解枚举的 i,j,k 
			{
				comp u=a[k],w={cos(o*Pi*K/i),sin(o*Pi*K/i)},v=a[k+i]*w;
				a[k]=u+v,a[k+i]=u-v;//转移文中写的 i,i+n/2 
			}
		}
	}
}
int main()
{
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i].x;
	for(int i=0;i<=m;i++) cin>>a[i].y;mmax=bger(n+m);
	DFT(a,mmax,1);
	for(int i=0;i<mmax;i++) a[i]=a[i]*a[i];
	DFT(a,mmax,-1);
	for(int i=0;i<=n+m;i++) cout<<(int)(a[i].y/2/mmax+0.5)<<" ";
	return 0;
}

同样优化显著。

NTT#

简单程度鉴定为：理解 $\texttt{FFT}$ 就能理解的，代码对着 $\texttt{FFT}$ 改即可。

快速计算 $a\times b=c$，其中 $c_k=\sum\limits_{i+j=k} a_ib_j$。结果对 $998244353$ 取模。模板题依然可以用P3803。

注意到 $\texttt{FFT}$ 中 $w_n$ 可以用模素数下的原根代替，其中由于补齐到 $2$ 的幂次了，$n$ 可以表示为 $2^k$。

具体的：设 $g$ 为模素数 $p$ 下的原根，其中 $p$ 为题目给定的模数。此时取 $w_{2^k}=g^{(p-1)/2^k}$ 则满足上文单位根需要的所有性质。

注意：此时素模数 $p$ 要满足，对于最高的 $2^k$，有 $2^k\mid p-1$。

常见模数 $998244353=119\times 2^{23}+1$，有原根 $3$。$167772161=5\times 2^{25}+1$，有原根 $3$。原根表 1，原根表 2。

直接把所有复数运算替换成简单的加法乘法取模即可。然后原根幂次要处理一下。常数和代码量均减少许多。

尽管此时要三次 DFT/IDFT，但由于不需要复数类速度依然吊打 FFT。

$\texttt{code}$

#include<bits/stdc++.h>
#define LL long long
//#define LD double
using namespace std;
const int N=4e6+5,mod=998244353,g=3,ig=332748118;//ig 为 3 的逆元 
//const LD Pi=acos(-1.0);
int mmax,n,m,rev[N],a[N],b[N];
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}//bger(n) 表示 >n 的第一个 2 的幂次
void DFT(int *a,int n,int o)
{
	for(int i=1;i<n;i++) rev[i]=rev[i>>1]>>1,(i&1)&&(rev[i]|=n/2);//求 rev 
	for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);//交换，注意别交换两次了 
	for(int i=1;i<n;i<<=1)
	{
		int s=ksm(o?g:ig,(mod-1)/(i<<1));//表示 w_{2*i} 单位根所对应的数 
		for(int j=0;j<n;j+=(i<<1))
		{
			for(int k=j,w=1;k<j+i;k++,w=1ll*w*s%mod)//w 循环遍历了 w_{2*i}^K
			{
				int u=a[k],v=1ll*a[k+i]*w%mod;
				a[k]=md(u+v),a[k+i]=md(u-v+mod);
			}
		}
	}
}
int main()
{
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i];
	for(int i=0;i<=m;i++) cin>>b[i];mmax=bger(n+m);
	DFT(a,mmax,1);DFT(b,mmax,1);
	for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
	DFT(a,mmax,0);int I=ksm(mmax,mod-2);
	for(int i=0;i<=n+m;i++) cout<<1ll*a[i]*I%mod<<" ";
	return 0;
}

为了更简洁（方便下面优化），对 IDFT 的算法进行一些小修改。

考虑 $(\star)$。有：$A^2:A^2_{i,j}=\begin{cases}n(n\mid i+j)\\ 0(n\nmid i+j)\end{cases}$。即 $A^2=\begin{bmatrix} n&0&0&\cdots &0& 0\\ 0&0&0&\cdots &0&n\\ 0&0&0&\cdots &n&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&0&n&\cdots &0& 0\\ 0&n&0&\cdots &0& 0\\ \end{bmatrix}$。$(\star\star)$

于是对 $a$ 做 IDFT 相当于先对 $a$ 做一次 DFT，翻转 $a_{1,2,\cdots ,n-1}$，然后对 $a_{0,2,\cdots ,n-1}$ 都 $/n$ 即可。这时就不需要 $g$ 的逆元了。

$\texttt{code}$

#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=4e6+5,mod=998244353,g=3;
int mmax,n,m,rev[N],a[N],b[N];
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}//bger(n) 表示 >n 的第一个 2 的幂次
void DFT(int *a,int n,int o)
{
	for(int i=1;i<n;i++) rev[i]=rev[i>>1]>>1,(i&1)&&(rev[i]|=n/2);//求 rev 
	for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);//交换，注意别交换两次了 
	for(int i=1;i<n;i<<=1)
	{
		int s=ksm(g,(mod-1)/(i<<1));//都用 g 来做 DFT 
		for(int j=0;j<n;j+=(i<<1))
		{
			for(int k=j,w=1;k<j+i;k++,w=1ll*w*s%mod)//w 循环遍历了 w_{2*i}^K
			{
				int u=a[k],v=1ll*a[k+i]*w%mod;
				a[k]=md(u+v),a[k+i]=md(u-v+mod);
			}
		}
	}
	if(!o)
	{
		reverse(a+1,a+n);int I=ksm(n,mod-2);
		for(int i=0;i<n;i++) a[i]=1ll*a[i]*I%mod;
	}//做 IDFT 要 rev 一下然后 /n 
}
int main()
{
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n>>m;
	for(int i=0;i<=n;i++) cin>>a[i];
	for(int i=0;i<=m;i++) cin>>b[i];mmax=bger(n+m);
	DFT(a,mmax,1);DFT(b,mmax,1);
	for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;
	DFT(a,mmax,0);
	for(int i=0;i<=n+m;i++) cout<<a[i]<<" ";
	return 0;
}

DFT/IDFT 的线性性#

注意到 DFT/IDFT 是线性变换，对于矩阵乘法形式，有：$A\vec{u}+A\vec{v}=A(\vec{u}+\vec{v})$。

考虑一种情况：要求 $f(x)g(x)+F(x)G(x)$ 的各项系数。若直接两次多项式乘法相加需要 $4$ 次 DFT，两次 IDFT。

而考虑 $4$ 次 DFT 后变成了 $u,v,U,V$，此时令 $u_i=u_iv_i+U_iV_i$，其中 $u_i$ 表示 $[x^i]u(x)$。然后对 $u$ 做一次 IDFT 即可。

此时减少了一次 IDFT！

---

考虑若求 $f(x)g(x)h(x)$，则可以对三个多项式做一次 DFT，然后把三处点值相乘，再做一次 IDFT 即可。

通过线性性可以减少许多冗余的 DFT/IDFT。

Super NTT(DIT-DIF)#

超级快速科技啊！科技目的为省去蝴蝶变换，优化常数。能写出飞速的 NTT，可以再也别担心多项式题卡常了！

参考文章 1，参考文章 2，有些大神看这里的转移图就能理解了。

下文记 $H(a)$ 为蝴蝶变换。显然这是线性变换，写成矩阵形式为：$\vec{u}\to H\vec{u}$，其中 $H$ 矩阵满足 $H_{i,rev(i)}=1$。

注意由于蝴蝶变换的性质，有：$H^{-1}=H,H^2=I$，其中 $I$ 为单位矩阵。

考虑从 IDFT 开刀。若令 $\texttt{IDFT}'(a)=\texttt{IDFT}(H(a))$，则按照先前 IDFT 流程做就可以免去蝴蝶变换。

记此时 DFT 矩阵为 $C$，则 $C^{-1}=A^{-1}\times H$，于是 $C=HA$。

则 $C_{i,j}=\sum\limits_{k=0}^{n-1} H_{i,k}A_{k,j}=A_{rev(i),j}$，第二个等号是由于只有 $H_{i,rev(i)}=1$。

验证下 $C$ 是否满足 DFT 矩阵条件：$C_{i,j+k}=A_{rev(i),j+k}=A_{rev(i),j}\times A_{rev(i),k}=C_{i,j}\times C_{i,k}$，显然是满足的。

考虑如何快速计算此时的 $\texttt{DFT}'(a)$，$\texttt{IDFT}'(a)$ 请自行推导。

$(C\vec{u})_i=\sum\limits_{j=0}^{n-1} C_{i,j}u_j=\sum\limits_{j=0}^{n-1} w_n^{rev(i)j}u_j=\sum\limits_{j=0}^{n/2-1} (w_n^{rev(i)j}u_j+w_n^{rev(i)(j+n/2)}u_{j+n/2})=\sum\limits_{j=0}^{n/2-1} w_n^{rev(i)j}(u_j+(-1)^{rev(i)}u_{j+n/2})$。

当 $0\le i<n/2$ 时，$2\mid rev(i),2\nmid rev(i+n/2)$。有：

$(C\vec{u})_i=\sum\limits_{j=0}^{n/2-1} w_n^{rev(i)j}(u_j+u_{j+n/2})$

$(C\vec{u})_{i+n/2}=\sum\limits_{j=0}^{n/2-1} w_n^{rev(i+n/2)j}(u_j-u_{j+n/2})=\sum\limits_{j=0}^{n/2-1} w_n^{(rev(i)+1)j}(u_j-u_{j+n/2})=\sum\limits_{j=0}^{n/2-1} w_n^{rev(i)j}\times w_n^{j}(u_j-u_{j+n/2})$

于是你从 $n\to n/2\to n/4$ 这样枚举，然后枚举对应的 $i,i+n/2$ 进行 $\left(u_i,u_{i+n/2}\right)\to\left(u_i+u_{i+n/2},w_n^{i}(u_i-u_{i+n/2})\right)$ 变换，然后就能递归到下一层做了。也省去了蝴蝶变换！

把变换 $C\vec{u}$ 称为做 DIF-FFT，变换 $C^{-1}\vec{u}$ 称为 DIT-FFT，此时省去了两次蝴蝶变换。同时原来 DFT 变换的所有性质也是几乎保留的。

---

考虑在枚举单位根的时候有 $O(n\log n)$ 次乘法取模，考虑优化。由于上文修改了 IDFT 形式，我们只需处理 $g$ 的若干幂次即可，不需要考虑 $g^{-1}$。

注意到有用的单位根形如：$w_n^{k}$，其中 $n$ 为 $2$ 的幂次，$0\le k<n$。

于是我们把所有合法的 $w_n^k$ 放在 $n+k$ 位置，然后找单位根直接顺序遍历一个区间即可。

这样只需要 $O(n)$ 预处理加上 $O(n\log n)$ 次连续内存的访问下标，有优秀的常数。

$\texttt{code}$

//洛谷 P3803
//https://www.luogu.com.cn/problem/P3803
#include<bits/stdc++.h>
#define fr(x) freopen(#x".in","r",stdin);freopen(#x".out","w",stdout);
using namespace std;
const int mod=998244353,N=4e6+5;
int n,m,a[N],b[N],w[N],mmax;
inline int rd()
{
    int x=0,zf=1;
    char ch=getchar();
    while(ch<'0'||ch>'9') (ch=='-')and(zf=-1),ch=getchar();
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
    return x*zf;
}
inline void wr(int x)
{
    if(x==0) return putchar('0'),putchar(' '),void();
    int num[35],len=0;
    while(x) num[++len]=x%10,x/=10;
    for(int i=len;i>=1;i--) putchar(num[i]+'0');
    putchar(' ');
}
inline int bger(int x){return x|=x>>1,x|=x>>2,x|=x>>4,x|=x>>8,x|=x>>16,x+1;}
inline int md(int x){return x>=mod?x-mod:x;}
inline int ksm(int x,int p){int s=1;for(;p;(p&1)&&(s=1ll*s*x%mod),x=1ll*x*x%mod,p>>=1);return s;}
//和之前一样的三个优化常数的函数
inline void init(int mmax)
{
	for(int i=1,j,k;i<mmax;i<<=1)
		for(w[j=i]=1,k=ksm(3,(mod-1)/(i<<1)),j++;j<(i<<1);j++)
			w[j]=1ll*w[j-1]*k%mod;
}//按照上面所述的方法预处理单位根
inline void DNT(int *a,int mmax)
{
	for(int i,j,k=mmax>>1,L,*W,*x,*y,z;k;k>>=1)//n->n/2->... 的枚举顺序
		for(L=k<<1,i=0;i<mmax;i+=L)
			for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
			//用了指针优化访问内存，注意 k 枚举的是 n/2，各个指针代表的含义请自行理解
				*y=1ll*(*x+mod-(z=*y))* *W%mod,*x=md(*x+z);
}//做 DFT 矩阵为 C 的 DFT 变换
inline void IDNT(int *a,int mmax)
{
	for(int i,j,k=1,L,*W,*x,*y,z;k<mmax;k<<=1)
		for(L=k<<1,i=0;i<mmax;i+=L)
			for(j=0,W=w+k,x=a+i,y=x+k;j<k;j++,W++,x++,y++)
				z=1ll* *W* *y%mod,*y=md(*x+mod-z),*x=md(*x+z);
	//和之前 IDFT 没啥区别，就是少了个蝴蝶变换
	//各个变量的含义和 DFT 时区别不大，区别在于 i,i+n/2 位置转移的系数
	reverse(a+1,a+mmax);//记得 rev
	for(int inv=ksm(mmax,mod-2),i=0;i<mmax;i++) a[i]=1ll*a[i]*inv%mod;//记得 /n
}//做 IDFT 矩阵为 C^{-1} 的 IDFT 变换
inline void NTT(int *a,int *b,int n,int m)
{
	mmax=bger(n+m);init(mmax);DNT(a,mmax);DNT(b,mmax);
	for(int i=0;i<mmax;i++) a[i]=1ll*a[i]*b[i]%mod;IDNT(a,mmax);
}//封装
int main()
{
	n=rd();m=rd();
	for(int i=0;i<=n;i++) a[i]=rd();
	for(int i=0;i<=m;i++) b[i]=rd();
	NTT(a,b,n,m);
	for(int i=0;i<=n+m;i++) wr(a[i]);
    return 0;
}

跑得飞快，足够应对所有多项式题了。

练习#

在此博客中有一大车总结的练习。拉一道出来讲一下：

P7844，求长为 $2^k$ 的数组在 $k$ 次 DFT 变换后的结果。即求 $A^k\vec{u}$，其中 $A$ 为范德蒙德矩阵。

考虑 $(\star\star)$。有：$A^2=\begin{bmatrix} n&0&0&\cdots &0& 0\\ 0&0&0&\cdots &0&n\\ 0&0&0&\cdots &n&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&0&n&\cdots &0& 0\\ 0&n&0&\cdots &0& 0\\ \end{bmatrix}$。

于是有：$A^4=(A^2)^2=n^2I_n$，其中 $I_n$ 为单位矩阵。

于是 $A^k\vec{u}=(n^2)^{\lfloor k/4\rfloor}A^{k\bmod 4}\vec{u}$，于是复杂度同 DFT 的复杂度，为 $O(n\log n)$。

FWT#

快速计算 $a\circ b=c$，其中 $c_k=\sum\limits_{i\circ j=k} a_ib_j$，需要做 $\circ=\texttt{and},\texttt{or},\texttt{xor}$ 的这些位运算的情况。模板题。

由于位运算，下文一般设数组长度为 $2^n$，不足补齐即可。下文 $\subseteq$ 一般表示二进制下的属于。

与、或时的 FWT#

考虑一般性推导，当 $\circ=\texttt{or}$ 时，注意到 $[j\ \texttt{or}\ k\subseteq i]=[j\subseteq i][k\subseteq i]$，于是取 $A_{i,j}=[j\subseteq i]$ 即可。

同理，当 $\circ=\texttt{and}$ 时，取 $A_{i,j}=[i\subseteq j]$ 即可。

---

考虑快速做 DFT 变换，拿 $\texttt{or}$ 举例，考虑快速求 $A\vec{u}$，即要快速求出 $b:b_i=\sum\limits_{j\subseteq i} a_j$。

考虑 $0\sim n-1$ 枚举二进制位，设枚举到第 $k$ 位，所有第 $k$ 位为 $0$ 的值 $i$，看成做 dp 转移 $b'_{i+2^k}\gets b_{i+2^k}+b_i$，初值 $b=a$ 即可，复杂度 $O(n2^n)$。

注：不明白可以想象 $b_{k,i}$ 表示对于前 $k$ 位 $j\subseteq i$ 的所有 $a_j$ 的和，然后做 dp 转移，发现可以简洁得做成上述过程。

inline void OR(modint *f) {
	for (int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
		for (int i = 0; i < n; i += o)
			for (int j = 0; j < k; j++)
				f[i+j+k] += f[i+j];
}

---

考虑快速做 IDFT 变换，即对上述数组知 $b$ 求 $a$，只需修改转移：$b'_{i+2^k}\gets b_{i+2^k} {\color{red}{\ -\ }} b_i$，初值依然不变，复杂度 $O(n2^n)$。

inline void IOR(modint *f) {
	for (int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
		for (int i = 0; i < n; i += o)
			for (int j = 0; j < k; j++)
				f[i+j+k] -= f[i+j];
}

异或时的 FWT#

此时不好直接构造 DFT 变换。由于位运算的按位独立性，尝试对 $i,j,k$ 二进制拆位，分别拆成 $p_{0,1,\cdots,n-1};q_{0,1,\cdots,n-1};r_{0,1,\cdots,n-1}$。

此时与、或情况下的 $A_{i,j}$ 可以表示为 $\prod\limits_{t=0}^{n-1} f_{p_i,q_i}$，其中 $f$ 为 $2\times 2$ 的矩阵。猜测 xor 时的 $A$ 矩阵也可以如此拆位表示！

观察此时 $f$ 需要满足的条件，即 $f_{p,q}f_{p,r}=f_{p,q\circ r}$，此时 $\circ=\texttt{xor}$。显然 $f$ 同样需要满足 $f$ 存在逆矩阵。

对 $q,r$ 取遍 $\{0,1\}^2$ 列出方程组：$\begin{cases} f_{p,0}^2=f_{p,0} \\ f_{p,1}^2=f_{p,0} \\ f_{p,0}f_{p,1}=f_{p,1} \end{cases}$

解得：$\begin{cases} f_{p,0}=0\\ f_{p,1}=0 \end{cases}$ 或 $\begin{cases} f_{p,0}=1\\ f_{p,1}=\pm 1 \end{cases}$

由于 $f$ 可逆，于是：前者可以舍去，而后者当 $p$ 分别取 $0,1$ 时 $f_{p,1}$ 的值应该不同。于是 $f=\begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}$ 或 $\begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}$。

由于前者关于对角线对称更加优美，于是选用前者。此时有 $f^{-1}=\dfrac{1}{2}f$。

注：可以思考当 $f=\begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix}$ 时如何模仿所述流程做，口胡一下。

此时可以刻画：$f_{p,q}=(-1)^{p\ \texttt{and}\ q}$，带入 $A$ 中：$A_{i,j}=(-1)^{\text{popcount}(i\ \texttt{and}\ j)}$，其中 $\text{popcount}$ 表示二进制位为 $1$ 的个数。

---

考虑快速做 DFT 变换，类似的，求 $A\vec{u}$，即要快速求出 $b:b_i=\sum\limits_{j}(-1)^{\text{popcount}(i\ \texttt{and}\ j)}a_j$。

设枚举到第 $k$ 位，转移：$\begin{cases} b'_{i}\gets b_i+b_{i+2^k}\\ b'_{i+2^k}\gets b_i-b_{i+2^k}\end{cases}$，即考虑 $\text{popcount}$ 的变化即可。IDFT 转移时乘个 $\dfrac{1}{2}$ 系数即可。

inline void XOR(modint *f, modint x = 1) {// x 取 1 和 1/2 时分别表示 DFT 和 IDFT 变换
	for (int o = 2, k = 1; o <= n; o <<= 1, k <<= 1)
		for (int i = 0; i < n; i += o)
			for (int j = 0; j < k; j++)
				f[i+j] += f[i+j+k],
				f[i+j+k] = f[i+j] - f[i+j+k] - f[i+j+k],
				f[i+j] *= x, f[i+j+k] *= x;
}

在 $\circ$ 为位运算时，一般也称 $f$ 为 DFT 矩阵。

注：注意到上文中的 $f$ 矩阵只是为了刻画 $A$ 优美形式的过渡式子，但看着这样简单的 $2\times 2$ 矩阵形式就应该有更大的用处。详见「其他进制下的其他位运算」。

OR/AND/XOR 卷积练习#

类似的，位运算的 DFT（或者称为 FWT） 也具有线性性，部分练习需要用到。

子集卷积，请自行查看题解。

一些直接的应用：P3175，CF1208F，CF1034E，CF1713F，ARC100C，CF1234F。

妙妙题。AGC034F，需要用到 FWT 的线性性推式子。

更多练习可以在学习子集卷积、集合幂级数后搜索。参考 cmd 博客的最后一个部分。

其他进制下的其他位运算#

再次回顾下 $\texttt{xor}$ 时的 dp 转移，发现即为：$\begin{bmatrix} b'_{i} \\ b'_{i+2^k} \end{bmatrix}=f\begin{bmatrix} b_{i} \\ b_{i+2^k} \end{bmatrix}$。而 $\texttt{and,or}$ 的情况也是如此。

考虑设 $2^k$ 时的 DFT 矩阵 $A$ 为 $A_k$。有 $f=A_1$，$A_k=\begin{bmatrix} f_{0,0}A_{k-1} & f_{0,1}A_{k-1} \\ f_{1,0}A_{k-1} & f_{1,1}A_{k-1} \end{bmatrix}$。

而 $\begin{bmatrix} b_{i} \\ b_{i+2^k} \end{bmatrix}$ 相当于在 $\times A_{k-1}$ 中提取出来，而 $\begin{bmatrix} b'_{i} \\ b'_{i+2^k} \end{bmatrix}$ 相当于在 $\times A_{k}$ 中提取出来。

于是 $\begin{bmatrix} b'_{i} \\ b'_{i+2^k} \end{bmatrix}=f\begin{bmatrix} b_{i} \\ b_{i+2^k} \end{bmatrix}$ 就容易理解了。

同时这也说明：对于一般进制的一般位运算，只需找到可逆矩阵 $f=A_1$，就可以用类似 $\texttt{or,and,xor}$ 的方法做转移，完成 DFT/IDFT。其中 IDFT 的转移矩阵自然就是 $f^{-1}$。

一些练习#

ABC288G，相当于要对 $a$ 做特殊位运算下的 IDFT，而此题 $f$ 矩阵时容易找到的，求下逆矩阵按上述流程做即可。

$\texttt{code}$

const int N=6e5+5;
int n,pw[N],a[N],b[3][3],A[3];
inline void FWT(int *A)
{
	static int B[3];B[0]=B[1]=B[2]=0;
	for(int i=0;i<3;i++) for(int j=0;j<3;j++) B[i]+=A[j]*b[i][j];
	for(int i=0;i<3;i++) A[i]=B[i];
}//矩阵成向量，按 IDFT 流程转移
int main()
{
	ios::sync_with_stdio(0);cin.tie(0);cout.tie(0);cin>>n;
	for(int i=pw[0]=1;i<=n;i++) pw[i]=3*pw[i-1];
	for(int i=0;i<pw[n];i++) cin>>a[i];
	b[0][0]=0;b[0][1]=1;b[0][2]=-1;
	b[1][0]=1;b[1][1]=-1;b[1][2]=1;
	b[2][0]=-1;b[2][1]=1;b[2][2]=0;//b 即为求得的 f 的逆矩阵
	for(int i=0;i<n;i++) for(int j=0;j<pw[n];j++) if(!((j/pw[i])%3))
	{
		for(int k=0;k<3;k++) A[k]=a[j+k*pw[i]];FWT(A);
		for(int k=0;k<3;k++) a[j+k*pw[i]]=A[k];
	}
	for(int i=0;i<pw[n];i++) cout<<a[i]<<" ";
	return 0;
}

gym102129A，思考当找不到（或不好找）可逆的 $f$ 时有啥处理方法，详见题解。

k 进制异或(不进位加法)#

考虑当 $\circ$ 表示 $k$ 进制异或（不进位加法）时如何求卷积。

考虑 $k=3$ 时的 $f$，求得为：$\begin{gather*} \begin{bmatrix} 1 & 1 & 1 \\ 1 & \omega & \omega^2 \\ 1 & \omega^2 & \omega^4 \\ \end{bmatrix} \end{gather*}$，其中 $\omega=w_3$，即 $3$ 次单位根。

发现这即是 $k$ 时的范德蒙德矩阵，证明按 $f_{p,q}f_{p,r}=f_{p,q\circ r}$ 直接代数计算即可。

对长度为 $k^n$ 的数组，复杂度为 $O(nk^n\times k^2)$，其中 $n$ 为枚举位数，$k^2$ 为矩阵乘向量的复杂度。

---

下文为作者还不会的，todo。

注意有时单位根 $w_k$ 在模意义下不存在，考虑造分圆多项式为“单位根”，参考 cmd 博客中「每一维加起来 mod k」部分。

例题：CF1103E Radix sum，P5577。