快速傅里叶变换
把多项式在系数表示与点值表示之间 (O(nlogn)) 转换
正变换是由系数 (Rightarrow) 点值
一个 (n) 次函数可用 (n+1) 个点值唯一确定,所以重点是选 (n+1) 个点快速求值
单位复数根
称 (x^n=1) 在复数意义下的解为 (n) 次复数根,显然有 (n) 个解。
联想高中学过的复数三角形式 (e^{ui}=cos(u)+isin(u)):
若 ((e^{ui})^n=1) ,则 (e^{uni}=cos(nu)+isin(nu)=1),则 (cos(nu)=1,sin(nu)=0 Rightarrow nu=2kpi,kin mathbb{Z})
于是 (u=frac{2kpi}{u},kin mathbb{Z}),所以 (e^{frac{2kpi i}{n}},kin mathbb{Z}) 都是解,当这无穷个数中其实只有 (n) 个不同的数(因为 (e^{ui}=e^{(u+2pi)i}))
设 (omega_{n}=e^{frac{2pi i}{n}}) ,则 (omega_n^k=e^{frac{2kpi i}{n}}),显然原方程的解是 ({omega_n^k|k=0,1,dots,n-1 })
称 (omega_n) 为单位复数根,下面是它的一些性质:
- (omega_n^n=1)
- 消去引理:(omega_{dn}^{dk}=omega_n^k)
证:显然 (omega_{dn}^{dk}=e^{frac{2dkpi i}{dn}}=e^{frac{2kpi i}{n}}=omega_n^k) - 折半引理:(n) 个 (n) 次单位复数根的平方的集合为 (n/2) 个 (n/2) 次单位复数根集合
证:由消去引理 (omega_n^{2k}=omega_{n/2}^k)
对于 (k=0,1,dots,n/2-1) ,恰好构成的是 ({omega_{n/2}^k| k=0,1,dots,n/2-1})
而对于 (k=k=n/2,n/2+1,dots,n-1),构成的是 ({omega_{n/2}^{k+n/2}| k=0,1,dots,n/2-1})
注意到 (omega_{n}^{k+n}=omega_n^nomega_n^k=omega_n^k) ,所以这两个集合是一样的 - (omega_{2n}^{n+k}=omega_{2n}^{k}omega_{2n}^{n}=-omega_{2n}^{k})
选择 ({omega_n^k|k=0,1,dots,n-1 }) ((n) 为2的幂) 为要计算点值的点
DFT
分治计算。
原多项式 (A(x)=sumlimits_{i=0}^{n-1} a_ix^i)
设 (A^{[0]}(x)=sumlimits_{i=0}^{n/2-1} a_{2i}x^{i},A^{[1]}(x)=sumlimits_{i=0}^{n/2-1} a_{2i+1}x^{i})
相当于把奇偶项分开了,则 (A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2))
由单位复数根的性质(公式中 (k<n/2) ) :
于是可以递归再合并了
当然代码里写的是非递归
IDFT
理论不太会……但代码与 (DFT) 几乎一样……
注意 (l) 的大小很玄妙……一定一定别开小了……
求 (r[]) :for(int i=1;i<l;i++) r[i]=(r[i>>1]>>1)|((i&1)*(l>>1));
(IDFT) 后的最终结果是 (int(A[i].r/l+0.5))
const double pi = 3.1415926535897932384626433832795;
struct c{ //复数
double r,i;
c() { r=i=0.0; }
c(double x,double y) { r=x; i=y; }
c operator + (const c &b) { return c(r+b.r,i+b.i); }
c operator += (const c &b) { return *this=*this+b; }
c operator - (const c &b) { return c(r-b.r,i-b.i); }
c operator -= (const c &b) { return *this=*this-b; }
c operator * (const c &b) { return c(r*b.r-i*b.i,r*b.i+i*b.r); }
c operator *= (const c &b) { return *this=*this*b; }
}x[N];
int l,r[N];
void fft(c *A,int ty){
for(int i=0;i<l;i++) x[r[i]]=A[i];
for(int i=0;i<l;i++) A[i]=x[i];
for(int i=2;i<=l;i<<=1){
c wn(cos(pi*2/i),ty*sin(pi*2/i));
for(int j=0;j<l;j+=i){
c w(1,0);
for(int k=j;k<j+i/2;k++){
c t=A[k+i/2]*w;
A[k+i/2]=A[k]-t;
A[k]+=t;
w*=wn;
}
}
}
}
快速数论变换
理论跟 (FFT) 差不多,区别是取模,将单位复数根换成了原根的幂
常见原根 (g) :998244353,1004535809的原根都是3
然而具体的我也不太清楚 (qwq)
int l,r[N],X[N];
void ntt(int *A,int ty){
for(int i=0;i<l;i++) X[r[i]]=A[i];
for(int i=0;i<l;i++) A[i]=X[i];
for(int i=2;i<=l;i<<=1){
int wn=Pow_mod(3,(P-1)/i);
if(ty==-1) wn=Pow_mod(wn,P-2);
for(int j=0;j<l;j+=i){
int w=1;
for(int k=j;k<j+i/2;k++){
int t=1ll*A[k+i/2]*w%P;
A[k+i/2]=Minus(A[k],t);
A[k]=Plus(A[k],t);
w=1ll*w*wn%P;
}
}
}
if(ty==1) return;
int Inv=Pow_mod(l,P-2);
for(int i=0;i<l;i++) A[i]=1ll*A[i]*Inv%P;
}
快速沃尔什变换
也是类似系数表示与点值表示之间的转换
与(&)卷积
(C[i]=sumlimits_{j&k=i} A[j] imes B[k])
需要设计出数组 (A'=FWT(A)) 满足 (C'[i]=A'[i] imes B'[i])
考虑这样设计 (A'[i]=sumlimits_{j&i=i} A[j])
那么 (C'[i]=sumlimits_{j&i=i} sumlimits_{k&l=j} A[k] imes B[l]=sumlimits_{j&i=i} A[j] imes sumlimits_{j&i=i} B[j]=A'[i] imes B'[i])
验证成功。
考虑如何求 (A'[]) (即 (FWT) 过程)
设 (A0[]) 表示下标最高位为0的,(A1[]) 表示下标最高位为1的
(FWT(A)=FWT(A0,A1)=(FWT(A0)+FWT(A1),FWT(A1)))
(IFWT(A)=IFWT(A0,A1)=(IFWT(A0)-IFWT(A1),IFWT(A1)))
void fwt_and(int *A,int ty){
for(int i=2;i<=n;i<<=1)
for(int j=0;j<n;j+=i)
for(int k=j;k<j+i/2;k++){
if(ty==1) A[k]=Plus(A[k],A[k+i/2]);
else A[k]=Minus(A[k],A[k+i/2]);
}
}
或(|)卷积
(C[i]=sumlimits_{j|k=i} A[j] imes B[k])
需要设计出数组 (A'=FWT(A)) 满足 (C'[i]=A'[i] imes B'[i])
考虑这样设计 (A'[i]=sumlimits_{j|i=i} A[j])
那么 (C'[i]=sumlimits_{j|i=i} sumlimits_{k|l=j} A[k] imes B[l]=sumlimits_{j|i=i} A[j] imes sumlimits_{j|i=i} B[j]=A'[i] imes B'[i])
验证成功。
考虑如何求 (A'[]) (即 (FWT) 过程)
设 (A0[]) 表示下标最高位为0的,(A1[]) 表示下标最高位为1的
(FWT(A)=FWT(A0,A1)=(FWT(A0),FWT(A0)+FWT(A1)))
(IFWT(A)=IFWT(A0,A1)=(IFWT(A0),IFWT(A1)-IFWT(A0)))
void fwt_or(int *A,int ty){
for(int i=2;i<=n;i<<=1)
for(int j=0;j<n;j+=i)
for(int k=j;k<j+i/2;k++){
if(ty==1) A[k+i/2]=Plus(A[k+i/2],A[k]);
else A[k+i/2]=Minus(A[k+i/2],A[k]);
}
}
异或(^)卷积
这个有点难……
原理:令 (i&j) 中 1 的个数的奇偶性为 “(i) 与 (j) 的奇偶性”,那么 “(i) 与 (k) 的奇偶性” 异或 “(k) 与 (j) 的奇偶性” 等于 “(i hat{} j) 与 (k) 的奇偶性”
令 (A'[i]=sumlimits_{i&j奇偶性为0} A[j]-sumlimits_{i&j奇偶性为1} A[j])
则
验证成功。
考虑如何求 (A'[]) (即 (FWT) 过程)
设 (A0[]) 表示下标最高位为0的,(A1[]) 表示下标最高位为1的
(FWT(A)=FWT(A0,A1)=(FWT(A0)+FWT(A1),FWT(A0)-FWT(A1)))
(IFWT(A)=IFWT(A0,A1)=(frac{IFWT(A0)+IFWT(A1)}{2},frac{IFWT(A0)-IFWT(A1)}{2}))
void fwt_xor(int *A,int ty){
for(int i=2;i<=n;i<<=1)
for(int j=0;j<n;j+=i)
for(int k=j;k<j+i/2;k++){
int t=A[k+i/2];
A[k+i/2]=Minus(A[k],t);
A[k]=Plus(A[k],t);
if(ty==-1)
A[k+i/2]=1ll*A[k+i/2]*Q%P,A[k]=1ll*A[k]*Q%P;//Q=1/2(mod P)
}
}
子集卷积
(C[i]=sumlimits_{j|k=i,j&k=0} A[j] imes B[k])
如果没有 (j&k=0) 的限制,就是或卷积。
有了这个限制就要求在 (j|k=i) 的基础上满足 (|i|=|j|+|k|) ,绝对值表示1的个数
而 (|i|=|j|+|k|) 让我们想到了普通的卷积。
于是我们卷积套卷积,外层普通卷积,内层或卷积。
具体地,令 (f[i][j]=[|j|==i] imes A[j],g[i][j]=[|j|==i] imes B[j])
那么每个 (f[i],g[i]) 都可看成是多项式
令 (h[i]=sumlimits_{j+k=i} f[j] oplus g[k]) ,两个多项式 (aoplus b) 表示或卷积
最终 (C[i]=h[|i|][i])
下面的是洛谷模板题代码
#include<cstdio>
#include<iostream>
#include<algorithm>
#define P 1000000009
using namespace std;
int read(){
int x=0;
char ch=getchar();
while(!isdigit(ch)) ch=getchar();
while(isdigit(ch)) x=x*10+ch-'0',ch=getchar();
return x;
}
const int N = 1<<20;
int n,m,a[N],b[N];
int Plus(int x,int y) { return x+y<P?x+y:x+y-P; }
int Minus(int x,int y) { return x>=y?x-y:x-y+P; }
void fwt(int *A,int ty){
for(int i=2;i<=m;i<<=1)
for(int j=0;j<m;j+=i)
for(int k=j;k<j+i/2;k++){
if(ty==1) A[k+i/2]=Plus(A[k+i/2],A[k]);
else A[k+i/2]=Minus(A[k+i/2],A[k]);
}
}
int f[25][N],g[25][N],num[N];
int h[25][N];
int main()
{
n=read();
m=(1<<n);
for(int i=0;i<m;i++) a[i]=read();
for(int i=0;i<m;i++) b[i]=read();
for(int i=0;i<m;i++){
int t=0,x=i;
for(;x;x>>=1) t+=(x&1);
num[i]=t;
f[t][i]=a[i];
g[t][i]=b[i];
}
for(int i=0;i<=n;i++){ //pre fwt
fwt(f[i],1); fwt(g[i],1);
}
for(int i=0;i<=n;i++)
for(int j=0;j<=i;j++)
for(int k=0;k<m;k++)
h[i][k]=Plus(h[i][k],1ll*f[j][k]*g[i-j][k]%P); //don't fwt back
for(int i=0;i<=n;i++) fwt(h[i],-1);
for(int i=0;i<m;i++) printf("%d ",h[num[i]][i]);
return 0;
}
狄利克雷卷积
(h=f*g),即 (h[i]=sumlimits_{d|i} f[d] imes g[frac{i}{d}])
在杜教筛中会用到
主要是有两个重要的式子,这里先声明几个函数:
(mu(n)) 莫比乌斯函数,(phi(n)) 欧拉函数
(epsilon(n)=[n=1]) 元函数,(I(n)=1) 恒等函数,(id(n)=n) 单位函数
重要式子来了:
- (mu * I=epsilon) 这个式子确保了莫比乌斯反演的正确性
- (phi * I=id) 看着很美妙但不会证…