【算法】快速傅里叶变换(FFT)初探
参考（大部分证明摘自）：https://oi.men.ci/fft-notes/
【简介】
快速傅里叶变换（FFT）是一种可以在$O(nlogn)$时间内完成的离散傅里叶变换（DFT）算法，在OI中主要用于加速向量卷积/多项式乘法运算。
【前置技能】
【引入】
有两个多项式$A(x)$和$B(x)$，求$C(x)=A(x)*B(x)$。
$A(x)=\sum_{i=0}^{n-1}a_ix^i$
$B(x)=\sum_{i=0}^{n-1}b_ix^i$
$C(x)=A(x)*B(x)=\sum_{i=0}^{2n-2}(\sum_{j=0}^i a_j*b_{i-j})x^i$
用朴素的做法复杂度显然是$O(n^2)$的。
如果再推推？
$C(x)=A(x)*B(x)$
$=\sum_{i=0}^{2n-2}(\sum_{j=0}^i a_j*b_{i-j})x^i$
$=\sum_{i=0}^{2n-2}\sum_{j=0}^i(a_j*x^j)(b_{i-j}*x^{i-j})$
$=\sum_{j=0}^{2n-2}(a_j*x^j)\sum_{i=0}^{2n-2}(b_i*x^i)$
式子左边和右边分别是多项式$A(x),B(x)$的点值表示，求一个多项式的点值表示和通过点值表示求多项式都可以通过FFT和IDFT做到$O(nlogn)$，所以我们就可以在$O(nlogn)$的时间进行多项式乘法运算了。
【系数表示法】
设$A(x)$表示一个$n-1$次多项式，所有项的系数组成的$n$维向量$(a_0, a_1, a_2, …, a_{n-1})$唯一确定了一个多项式，这种表示法叫系数表示法。
$A(x)=\sum_{i=0}^{n-1}a_ix^i$
【点值表示法】
如果将一组互不相同的$n$个$x$代入多项式，得到$n$个值，组成的$n$维向量为$(x_0, x_1, x_2, …, x_{n-1}),(y_0, y_1, y_2, …, y_{n-1})$，唯一确定了一个多项式（证明在下方），这种表示法叫点值表示法。
$y_i=A(x_i)=\sum_{j=0}^{n-1} a_jx_i^j$
一个$n-1$次多项式在$n$个不同点的取值唯一确定一个多项式的证明：
若命题不成立，则存在两个不同的$n-1$次多项式$A(x),B(x)$满足$\forall i\subseteq [0,n-1]，A(x_i)=B(x_i)$。
令$C(x)=A(x)-B(x)$，则$C(x)$也是一个$n-1$次多项式，且满足$\forall i\subseteq [0,n-1]，C(x_i)=0$。
即$C(x)$有$n$个根，与代数基本定理（一个$n-1$次多项式在复数域上至多只有$n-1$个根）相矛盾，故原命题成立。
【卷积】
$(f*g)(n)$称为$f,g$的卷积，其离散的定义为：
$(f*g)(n)=\sum_{-\infty }^{\infty}f(\tau)g(n-\tau)$
特征为$\tau+(n-\tau)=n$
【复数、复平面】相关知识见：http://www.ruanyifeng.com/blog/2012/09/imaginary_number.html
【单位根】
在复平面上以原点为起点，单位圆的$n$等分点为终点，作$n$个向量，所得的幅角为正且最小的向量对应的复数$\omega _n$，称为$n$次单位根，则其余的$n-1$个向量对应的复数分别为$\omega _n^2,\omega _n^3…,\omega _n^n$，显然$\omega _n^0=\omega _n^n=1$。
因为单位根的幅角为$\frac{1}{n}$，所以有$\omega _n^k=cos k\frac{2\pi}{n}+i sin k\frac{2\pi}{n}$。
【单位根的性质】
①$\omega _{2n}^{2k}=\omega _n^k$
②$\omega _n^{k+\frac{n}{2}}=-\omega _n^k$
结合单位圆较易理解，证明略
【快速傅里叶变换(FFT)】
点值向量$(A(\omega_n^0),A(\omega_n^1),A(\omega_n^2),…,A(\omega_n^{n-1}))$被称为其系数向量$(a_0, a_1, a_2, …, a_{n-1})$的离散傅里叶变换。
把$A(x)$按奇偶次项分组。
令
$A_1(x)=(a_0+a_2x+a_4x^2+…+a_{n-2}x^{\frac{n}{2}-1})$
$A_2(x)=(a_1+a_3x+a_5x^2+…+a_{n-1}x^{\frac{n}{2}-1})$
$A(x)=A_1(x^2)+xA_2(x^2)$
若$k<\frac{n}{2}$，求$A(\omega _n^k)$：
$A(\omega _n^k)=A_1(\omega _n^{2k})+\omega _n^kA_2(\omega _n^{2k})$
$=A_1(\omega _{\frac{n}{2}}^k)+\omega _n^kA_2(\omega _{\frac{n}{2}}^k)$
求$A(\omega _n^{k+\frac{n}{2}})$：
$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega _n^{2k+n})+\omega _n^{k+\frac{n}{2}}A_2(\omega _n^{2k+n})$
$=A_1(\omega _n^{2k}\times \omega_n^n)-\omega_n^kA_2(\omega_n^{2k}\times \omega_n^n)$
$=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})$
$=A_1(\omega_{\frac{n}{2}}^k)-\omega_n^kA_2(\omega_{\frac{n}{2}}^k)$
根据上方的推导，我们显然可以使用一个分治的算法来求得所有的点值，效率$O(nlogn)$。
这就是Cooley-Tukey算法。
【傅里叶逆变换(IDFT)】
将点值表示的多项式转换为系数表示，可以使用傅里叶逆变换。
设$(y_0, y_1, y_2, …, y_{n-1})$为$(a_0, a_1, a_2, …, a_{n-1})$的离散傅里叶变换，设另一个向量$(c_0, c_1, c_2, …, c_{n-1})$满足
$c_k=\sum_{i=0}^{n-1}y_i(\omega _n^{-k})^i$
$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}(a_j(\omega _n^i)^j)(\omega_n^{-k})^i$
$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}(a_j(\omega _n^j)^i(\omega_n^{-k})^i)$
$=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omega_n^{j-k})^i$
$=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)$
设
$S(\omega_n^k)=1+\omega_n^k+(\omega_n^k)^2+…+(\omega_n^k)^{n-1}$
$k=0$时，显然$S(\omega_n^k)=n$
$k\neq 0$时，则等式两边乘上$\omega_n^k$
$\omega_n^kS(\omega_n^k)=\omega_n^k+(\omega_n^k)^2+(\omega_n^k)^3+…+(\omega_n^k)^{n}$
两式相减，得
$\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^n-1$
$S(\omega_n^k)=\frac{(\omega_n^k)^n-1}{\omega_n^k-1}$
分子为0，分母不为0，故
$S(\omega_n^k)=0$
考虑上式
$c_k=\sum_{j=0}^{n-1}a_j(\sum_{i=0}^{n-1}(\omega_n^{j-k})^i)$
$=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k})$
当$j=k$时，$S(\omega_n^{j-k})=n$，否则为0，故
$c_i=na_i$
$a_i=\frac{1}{n}c_i$
所以使用单位根的倒数作为$x$，做一次FFT，对得到的每个数除以$n$，即可得到点值表示对应的系数表示。
【实现】
【迭代FFT】
递归版FFT速度十分慢，所以一般使用迭代版的FFT。
观察规律可知FFT分治到终止时，某个编号的初始下标是终止下标二进制位翻转后的值，然后就可以模拟分治进行迭代FFT了。
【蝴蝶操作】
运用一个临时变量省去一个中转数组。
若此时$a(k)$为$A_1(\omega_{\frac{n}{2}}^k)$，$a(k+\frac{n}{2})$为$A_2(\omega_{\frac{n}{2}}^k)$，用$a(k)+\omega_n^k\times a(\frac{n}{2}+k)$更新$a(k)$，用$a(k)-\omega_n^k\times a(\frac{n}{2}+k)$更新$a(\frac{n}{2}+k)$时，用一个临时变量$t$进行以下操作。
$t=\omega_n^k\times a(\frac{n}{2}+k)$
$a(\frac{n}{2}+k)=a(k)-t$
$a(k)=a(k)+t$
【注意】最后求系数取整时为减少误差一般+0.5以四舍五入
【流程】
将数从初始位置换到终止位置。
模拟分治进行迭代FFT。
【例题】
裸的大数乘法，显然两个大数可以看成两个多项式
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> using namespace std; const int maxn=500010, inf=1e9; const double pi=acos(-1); struct cpx{double r, i; cpx(double _r=0, double _i=0):r(_r),i(_i){};}a[maxn], b[maxn], c[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);} int n, mx; int xs[maxn]; char s[maxn]; inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } } int main() { scanf("%d", &n); for(mx=1;mx<(n<<1);mx<<=1); scanf("%s", s+1); for(int i=0;i<n;i++) a[i]=s[n-i]-\'0\'; scanf("%s", s+1); for(int i=0;i<n;i++) b[i]=s[n-i]-\'0\'; fft(mx, a, 1); fft(mx, b, 1); for(int i=0;i<mx;i++) c[i]=a[i]*b[i]; fft(mx, c, -1); for(int i=0;i<mx;i++) xs[i]=(int)(c[i].r/mx+0.5); for(int i=0;i<mx;i++) if(xs[i]>=10) { xs[i+1]+=xs[i]/10, xs[i]%=10; if(i==mx-1) mx++; } mx--; while(!xs[mx] && mx>0) mx--; for(int i=mx;~i;i--) printf("%d", xs[i]); }
View Code
$c[k]=\sum_{i=k}^{n-1} a[i]*b[i-k]$不是卷积的形式，但化一化把$a$倒过来就是卷积的形式了。
$c[k]=ans[n-k-1]=\sum_{i=0}^{n-k-1}a[n-k-i-1]*b[i]$，所以FFT加速后输出$ans[n-1]$~$ans[0]$即可。
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> using namespace std; const int maxn=500010, inf=1e9; const double pi=acos(-1); struct cpx{double r, i; cpx(double _r=0, double _i=0):r(_r),i(_i){};}a[maxn], b[maxn], c[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);} int n, x, mx; int xs[maxn]; inline void read(int &k) { int f=1; k=0; char c=getchar(); while(c<\'0\' || c>\'9\') c==\'-\'&&(f=-1), c=getchar(); while(c<=\'9\' && c>=\'0\') k=k*10+c-\'0\', c=getchar(); k*=f; } inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } } int main() { read(n); for(mx=1;mx<(n<<1);mx<<=1); for(int i=0;i<n;i++) read(x), a[n-i-1]=cpx(x), read(x), b[i]=cpx(x); fft(mx, a, 1); fft(mx, b, 1); for(int i=0;i<mx;i++) c[i]=a[i]*b[i]; fft(mx, c, -1); for(int i=0;i<n;i++) xs[i]=(int)(c[i].r/mx+0.5); for(int i=n-1;~i;i--) printf("%d\n", xs[i]); }
View Code
$E_i=\sum_{j=0}^{i-1}\frac{q_j}{(i-j)^2}-\sum_{j=i+1}^{n-1}\frac{q_j}{(j-i)^2}$
$A_i=q_i$
$B_i=\frac{1}{i^2}$
$E_i=\sum_{j=0}^{i-1}A_j*B_{i-j}-\sum_{j=i+1}^{n-1}A_j*B_{j-i}$
注意$B_0$取不到，需要设为0，然后$j=i$就可以取了，左边直接FFT。
右半部分不是卷积的形式，需要转化：
$\sum_{j=i+1}^{n-1}A_j*B_{j-i}$
$\sum_{j=1}^{n-i-1}A_{i+j}*B_{j}$
$\sum_{j=1}^{n-i-1}A’_{n-j-i-1}*B_{j}=Ans2_{n-i-1}$
$E_i=Ans1_i-Ans2_{n-i-1}$
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> #define ll long long using namespace std; const int maxn=500010, inf=1e9; const double pi=acos(-1); struct cpx{double r, i; cpx(double _r=0, double _i=0):r(_r),i(_i){};} a1[maxn], a2[maxn], b[maxn], ans1[maxn], ans2[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);} int n, mx; double x; inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } } int main() { scanf("%d", &n); for(mx=1;mx<(n<<1);mx<<=1); for(int i=0;i<n;i++) scanf("%lf", &x), a1[i]=cpx(x), a2[n-i-1]=cpx(x); for(int i=1;i<n;i++) b[i]=cpx(1.0/i/i); fft(mx, a1, 1); fft(mx, a2, 1); fft(mx, b, 1); for(int i=0;i<mx;i++) ans1[i]=a1[i]*b[i]; for(int i=0;i<mx;i++) ans2[i]=a2[i]*b[i]; fft(mx, ans1, -1); fft(mx, ans2, -1); for(int i=0;i<n;i++) printf("%.3lf\n", ans1[i].r/mx-ans2[n-i-1].r/mx); }
View Code
好厉害的题，完全没看出是FFT
不连续回文子序列=回文子序列-连续回文子序列
后者显然可以用manacher求得
对于回文子序列，只要知道了一个位置左右两边有多少对关于此位置对称的相同的字符$f_i$，用$2^{f_i}-1$就能得到答案了。
求$f_i$我们同样可以使用manacher的思想，先在原串$S$中插入一下分隔符，如$aba$改为$\$\# a\# b\# a\# $得到新串$S\’$。
考虑两个相等的字符对哪个位置有贡献，显然这三个位置组成了一个等差数列，故$x+y=2mid$，因为$x$和$y$是字符，且字符在$S\’$中下标为偶数，所以有$\frac{x}{2}+\frac{y}{2}=mid$，而$\frac{x}{2}$在原串刚好是字符位置$+1$的位置。那么就能知道，对于$S_x=S_y$，他们对$x+y+2$有贡献。
对$a$和$b$的贡献分别算，算$a$时把字符是$a$的位置贡献$gx_i$设为1，$b$的位置的贡献$gx_i$设为0，算$b$贡献时则反之。
那么我们就能得到一个求$f_i$的公式：
可以发现分子是一个卷积的形式，直接上FFT就好了。
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> #define MOD(x) ((x)>=mod?(x)-mod:(x)) #define ll long long using namespace std; const int maxn=500010, inf=1e9, mod=1e9+7; const double pi=acos(-1); struct cpx{double r, i;cpx(double _r=0, double _i=0):r(_r), i(_i){};} a[maxn], b[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);} int n, N, ans, mid, mx, mxx; int p[maxn], mi[maxn], f[maxn]; char ch[maxn], s[maxn]; inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } if(f==-1) for(int i=0;i<N;i++) a[i].r/=N; } inline void getans(int ty) { for(int i=0;i<n;i++) a[i]=cpx(ch[i]-\'a\'==ty, 0); for(int i=n;i<mx;i++) a[i]=cpx(0, 0); fft(mx, a, 1); for(int i=0;i<mx;i++) b[i]=a[i]*a[i]; fft(mx, b, -1); for(int i=2;i<N;i++) f[i]+=(int)(b[i-2].r+0.5); } int main() { scanf("%s", ch); n=strlen(ch); for(mx=1;mx<(n<<1);mx<<=1); for(int i=1;i<=n;i++) s[i<<1]=ch[i-1], s[(i<<1)|1]=\'#\'; N=(n<<1)+2; s[1]=s[N]=\'#\'; ans=mid=mxx=0; for(int i=1;i<=N;i++) { if(mxx>i) p[i]=min(p[(mid<<1)-i], mxx-i); else p[i]=1; while(s[i-p[i]]==s[i+p[i]]) p[i]++; if(p[i]+i>mxx) mxx=p[i]+i, mid=i; ans=MOD(ans+mod-(p[i]>>1)); } mi[0]=1; for(int i=1;i<N;i++) mi[i]=1ll*(mi[i-1]<<1)%mod; getans(0); getans(1); for(int i=2;i<N;i++) ans=MOD(ans+mi[(f[i]+1)>>1]-1); printf("%d\n", ans); }
View Code
母函数+容斥用FFT优化…
对于每一个物品$i$，使$A[a_i]$++$,B[2*a_i]$++$,C[3*a_i]$++，即$A,B,C$三个多项式代表选$1$个，$2$个，$3$个同样物品的方案。
对于选$1$个物品的方案，显然就是$A$了。
对于选$2$个物品的方案，$A^2$会多算$(x_i, x_i)$的情况，于是需要减去$B$，故方案为$\frac{A^2-B}{2}$。
对于选$3$个物品的方案，$A^3$会多算$(x_i, x_i, x_i),(x_i, x_i, y_i),(x_i, y_i, x_i), (y_i, x_i, x_i)$的情况，后三者的方案数刚好是$A*B*3$，但是$A*B*3$还包括$3$个$(x_i, x_i, x_i)$，而我们只需要减去一个，所以要多加上$2*C$，故方案为$\frac{A^3-A*B*3+2*C}{6}$。
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> #define ll long long using namespace std; const int maxn=500010, inf=1e9; const double pi=acos(-1); struct cpx{double r, i; cpx(double _r=0, double _i=0):r(_r), i(_i){};} a[maxn], b[maxn], c[maxn], a2[maxn], a3[maxn], ab[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);} int n, mx, mxx, x; int ans[maxn]; inline void read(int &k) { int f=1; k=0; char c=getchar(); while(c<\'0\' || c>\'9\') c==\'-\'&&(f=-1), c=getchar(); while(c<=\'9\' && c>=\'0\') k=k*10+c-\'0\', c=getchar(); k*=f; } inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } if(f==-1) for(int i=0;i<N;i++) a[i].r/=N; } int main() { read(n); for(int i=0;i<n;i++) read(x), ans[x]++, a[x].r++, b[x<<1].r++, c[x*3].r++, mxx=max(mxx, 3*x); for(mx=1;mx<mxx;mx<<=1); fft(mx, a, 1); fft(mx, b, 1); fft(mx, c, 1); for(int i=0;i<mx;i++) a[i]=(a[i]*a[i]*a[i]-cpx(3)*a[i]*b[i]+cpx(2)*c[i])*cpx(1.0/6)+(a[i]*a[i]-b[i])*cpx(1.0/2)+a[i]; fft(mx, a, -1); for(int i=1;i<mx;i++) if((int)(a[i].r+0.5))printf("%d %d\n", i, (int)(a[i].r+0.5)); }
View Code
设旋转了$j$位，亮度增加了$c$，则差异值为$\sum_{i=0}^{n-1}(x_{i+j}-y_i+c)^2$
把式子拆开之后可以发现除了$-x_{i+j}*y_i$之外都是常数项和一次项，都可以直接预处理。
对$x$倍增一次，然后把$x$倒过来就可以对每种颜色求$-x_{i+j}*y_i$了。
$\sum_{i=0}^{n-j-1}x_{i+j}*y_i$
$\sum_{i=0}^{n-j-1}x_{n-i-j-1}*y_i$
然后枚举旋转次数和颜色即可。
233
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> #define ll long long using namespace std; const int maxn=500010; const ll inf=1e18; const double pi=acos(-1); struct cpx{double r, i;cpx(double _r=0, double _i=0):r(_r), i(_i){};} a[maxn], b[maxn], c[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.r*b.i+a.i*b.r);} int n, m, x, mx; ll sqrx, sumx, sqry, sumy; ll xy[maxn], ans; inline void read(int &k) { int f=1; k=0; char c=getchar(); while(c<\'0\' || c>\'9\') c==\'-\'&&(f=-1), c=getchar(); while(c<=\'9\' && c>=\'0\') k=k*10+c-\'0\', c=getchar(); k*=f; } inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } if(f==-1) for(int i=0;i<N;i++) a[i].r/=N; } int main() { read(n); read(m); for(int i=0;i<n;i++) read(x), sqrx+=1ll*x*x, sumx+=x, a[i]=a[i+n]=cpx(x); for(int i=0;i<n;i++) read(x), sqry+=1ll*x*x, sumy+=x, b[i]=cpx(x); n<<=1; for(mx=1;mx<(n<<1);mx<<=1); reverse(a, a+n); fft(mx, a, 1); fft(mx, b, 1); for(int i=0;i<mx;i++) c[i]=a[i]*b[i]; fft(mx, c, -1); for(int i=0;i<n;i++) xy[i]=(ll)(c[n-i-1].r+0.5); n>>=1; ans=inf; for(int i=0;i<n;i++) for(int j=-m;j<=m;j++) ans=min(ans, sqrx+sqry+1ll*j*j*n+2ll*(-xy[i]+sumx*j-sumy*j)); printf("%lld\n", ans); }
View Code
例7：codeforces 528D. Fuzzy Search
字符串匹配问题使用FFT的话，应该对于每一种字符分别统计。把母串对于这种字符能够匹配的位置$i$，令$A_i=1$，否则为$0$，对模式串该字符位置$j$，令$B_j=1$，否则为$0$。
那么求每个位置能否匹配，等价于求：
$\sum_{l=0}^{n-1}\sum_{i=0}^{l+m-1}A_{l+i}\cdot B_{i}$
$\sum_{l=0}^{n-1}\sum_{i=0}^{l+m-1}A_{l+i}\cdot B_{m-i-1}$
卷积形式直接上FFT即可。
虽然$i$超过$m-1$时候$B$会越界，但是FFT处理的时候实际上不会越界，完全可以放心搞…
对于这题预处理出母串每个位置能匹配哪些字符即可然后直接上这种做法最后判断每个位置能匹配的次数是否为模式串长度即可。
#include<iostream> #include<cstring> #include<cstdlib> #include<cstdio> #include<cmath> #include<algorithm> #define ll long long using namespace std; const int maxn=1000010, inf=1e9; const double pi=acos(-1); struct cpx{double r, i;cpx(double _r=0, double _i=0):r(_r), i(_i){};} a[maxn], b[maxn], c[maxn]; cpx operator + (cpx a, cpx b){return cpx(a.r+b.r, a.i+b.i);} cpx operator - (cpx a, cpx b){return cpx(a.r-b.r, a.i-b.i);} cpx operator * (cpx a, cpx b){return cpx(a.r*b.r-a.i*b.i, a.i*b.r+a.r*b.i);} int n, m, k, N, ANS; int ok[maxn][4], cnt[4], ans[maxn]; char s[maxn], t[maxn]; inline void read(int &k) { int f=1; k=0; char c=getchar(); while(c<\'0\' || c>\'9\') c==\'-\'&&(f=-1), c=getchar(); while(c<=\'9\' && c>=\'0\') k=k*10+c-\'0\', c=getchar(); k*=f; } inline int col(char ch) { if(ch==\'A\') return 0; if(ch==\'T\') return 1; if(ch==\'G\') return 2; return 3; } inline void fft(int N, cpx a[], int f) { for(int i=0, j=0;i<N;i++) { if(i<j) swap(a[i], a[j]); for(int k=N>>1;(j^=k)<k;k>>=1); } for(int i=2;i<=N;i<<=1) for(int j=0, m=i>>1;j<N;j+=i) { cpx nw=cpx(cos(2.*pi/i), sin(2.*pi/i)*f), w=cpx(1, 0); for(int k=0;k<m;k++) { cpx t=a[j+k+m]*w; a[j+k+m]=a[j+k]-t; a[j+k]=a[j+k]+t; w=w*nw; } } if(f==-1) for(int i=0;i<N;i++) a[i].r/=N; } int main() { read(n); read(m); read(k); scanf("%s", s); scanf("%s", t); for(int i=0;i<n;i++) { cnt[col(s[i])]++; if(i-k-1>=0) cnt[col(s[i-k-1])]--; for(int j=0;j<4;j++) if(cnt[j]) ok[i][j]=1; } cnt[0]=cnt[1]=cnt[2]=cnt[3]=0; for(int i=n-1;~i;i--) { cnt[col(s[i])]++; if(i+k+1<n) cnt[col(s[i+k+1])]--; for(int j=0;j<4;j++) if(cnt[j]) ok[i][j]=1; } for(N=1;N<(n+m);N<<=1); for(int i=0;i<4;i++) { memset(a, 0, sizeof(a)); memset(b, 0, sizeof(b)); for(int j=0;j<n;j++) a[j]=cpx(ok[j][i], 0); for(int j=0;j<m;j++) b[m-j-1]=cpx(col(t[j])==i, 0); fft(N, a, 1); fft(N, b, 1); for(int j=0;j<N;j++) c[j]=a[j]*b[j]; fft(N, c, -1); for(int j=0;j<n;j++) ans[j]+=(int)(c[j+m-1].r+0.5); } for(int i=0;i<n;i++) if(ans[i]==m) ANS++; printf("%d\n", ANS); }
View Code
例8：