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)\)。
\((C\vec{u})_i=\sum\limits_{j=0}^{n-1} C_{i,j}u_j=\sum\limits_{j=0}^{n-1} w^{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\) 进行变换即可，也省去了蝴蝶变换！
把变换 \(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。
更多练习可以在学习子集卷积、集合幂级数后搜索。参考 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」部分。
