OI - TN's Kingdom III - Assassination(POJ)
作者：互联网
这个标题就出奇的长
# 题面
\(n\) 次多项式 \(\alpha\) 和 \(\beta\) 做循环卷积得到 \(\gamma\)。现在已知 \(\beta,\gamma\)，求 \(\alpha\)。
换句话说：
\[\gamma=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\alpha_i\beta_j\cdot x^{(i+j)\bmod n} \]数据规模：\(n< 2^{17}\)。
# 解析
直到现在我才知道 FFT 本来就是做的循环卷积……只不过以前用 FFT 做多项式线性卷积都是把多项式次数 \(N\) 设得特别大，然后体现在系数上就没有循环。
但是我们实现的 FFT 都是分治做的，要求多项式次数是二的幂。就没有办法做任意长度的循环卷积。
Hint.
普通的 FFT 也可以解决一部分循环卷积问题，只需要把多项式次数设得很大，先计算出 $\alpha,\beta$ 的线性卷积，再手动把 $x^t$ 的系数加到 $x^{t\bmod n}$ 上去。
只不过这样没有办法倒过来，知道 $\alpha\times\beta=\gamma$ 的 $\beta,\gamma$ 求 $\alpha$。
那就先不管 FFT，来推一下对多项式 \(f\) 求 DFT 的式子。记 \(X_k\) 表示 \(f(e^{-\frac{2\pi k}Ni})\)（也就是求 DFT 后的第 \(k\) 项），\(f_k\) 是 \(f\) 的 \(x^k\) 项的系数：
\[\begin{aligned} X_k&=\sum_{n=0}^{N-1}f_ne^{-\frac{2\pi kn}{N}i} \end{aligned} \]Hint.
$e^{\frac{2\pi} ni}$ 就是单位复根，也可以写作 $\omega_n=\cos\frac{2\pi}n+\sin\frac{2\pi}ni$（是不是看起来熟悉一点了）。
接下来就是最重要的部分：
\[\begin{aligned} &nk=\frac{(n-k)^2-n^2-k^2}2\\ &e^{-\frac{2\pi nk}Ni}=e^{\frac{-k^2\pi}N}\cdot e^{\frac{(k-n)^2\pi}{N}}\cdot e^{\frac{-n^2\pi}{N}} \end{aligned} \]于是
\[X_k=e^{-\frac{k^2\pi}{N}}\sum_{n=0}^{N-1}e^{\frac{(k-n)^2\pi}{N}}\cdot\big(e^{-\frac{n^2\pi}{N}}f_n\big) \]可以看出 \(X_k\) 后面的求和式是一个线性卷积，可以用 FFT。
值得注意的是这个卷积的下标 \(k-n\) 的范围是 \([1-N,N-1]\)，是可能有负数的。所以先给它平移 \(N\) 位。
逆变换也基本一样：
\[X_k'=e^{\frac{k^2\pi}{N}}\sum_{n=0}^{N-1}e^{-\frac{(k-n)^2\pi}{N}}\cdot\big(e^{\frac{n^2\pi}{N}}f_n\big) \]所以构造出这样两个数列：
\[A=\sum_{i=0}^{2N-1}e^{\frac{(i-N)^2\pi}{N}}x^i\\ B=\sum_{i=0}^{N-1}f_ie^{-\frac{i^2\pi}{N}}x^i \]卷积得到多项式 \(C\)，别忘了 \(C\) 是左移了 \(N\) 位的。
\[X_k=e^{-\frac{k^2\pi}{N}}C_{k+N} \]这样就实现了循环卷积。整个算法思路归为 bluestein，是一种把循环卷积转化为线性卷积的方法。
# 源代码
/*Lucky_Glass*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
struct Complex{
double r,i;
Complex(const double &_r=0,const double &_i=0):r(_r),i(_i){}
Complex operator *(const Complex &v)const{
return Complex(r*v.r-i*v.i,r*v.i+i*v.r);
}
Complex operator +(const Complex &v)const{
return Complex(r+v.r,i+v.i);
}
Complex operator -(const Complex &v)const{
return Complex(r-v.r,i-v.i);
}
Complex operator /(const double &k)const{
return Complex(r/k,i/k);
}
Complex operator /(const Complex &v)const{
return Complex((i*v.i+r*v.r)/(v.r*v.r+v.i*v.i),(i*v.r-r*v.i)/(v.r*v.r+v.i*v.i));
}
};
Complex omega(const double &k){return Complex(cos(k),sin(k));}
const int N=(2<<17)+10;
const double PI=acos(-1.0);
int rev[N<<2];
void FFT(Complex *ary,int n,int typ){
for(int i=1;i<n;i++){
rev[i]=rev[i>>1]>>1|((i&1)*(n>>1));
if(rev[i]<i) swap(ary[rev[i]],ary[i]);
}
for(int i=1,ii=2;i<n;i<<=1,ii<<=1){
Complex wn=omega(PI/i*typ);
for(int j=0;j<n;j+=ii){
Complex wi(1,0);
for(int k=j;k<j+i;k++,wi=wi*wn){
Complex tmp=wi*ary[k+i];
ary[k+i]=ary[k]-tmp;
ary[k]=ary[k]+tmp;
}
}
}
if(typ==-1) for(int i=0;i<n;i++) ary[i]=ary[i]/n;
}
void bluestein(Complex *ary,int n,int typ){
static Complex aryA[N<<2],aryB[N<<2];
memset(aryA,0,sizeof aryA);
memset(aryB,0,sizeof aryB);
for(int i=0;i<n;i++)
aryA[i]=omega(-typ*PI*i*i/n)*ary[i];
for(int i=0;i<(n<<1);i++)
aryB[i]=omega(typ*PI*(i-n)*(i-n)/n);
int len=1;
while(len<(n<<2)) len<<=1;
FFT(aryA,len,1),FFT(aryB,len,1);
for(int i=0;i<len;i++) aryA[i]=aryA[i]*aryB[i];
FFT(aryA,len,-1);
for(int i=0;i<n;i++){
ary[i]=aryA[i+n]*omega(-typ*PI*i*i/n);
if(typ==-1) ary[i]=ary[i]/n;
}
}
int n;
Complex aryA[N],aryB[N];
int main(){
// freopen("input.in","r",stdin);
scanf("%d",&n);
for(int i=0;i<n;i++) scanf("%lf",&aryA[i].r);
for(int i=0;i<n;i++) scanf("%lf",&aryB[i].r);
bluestein(aryA,n,1),bluestein(aryB,n,1);
for(int i=0;i<n;i++) aryA[i]=aryB[i]/aryA[i];
bluestein(aryA,n,-1);
for(int i=0;i<n;i++) printf("%.4f\n",aryA[i].r);
return 0;
}
THE END
Thanks for reading!
\[\begin{split} “\ &你是沙鸥\ 你是滴漏\\ &你是白昼\ 你是不朽\\ &你是从来都没有\ 缄默的口\\ &你是跌落深谷里\ 最后的柔\\ &你是行走于人间的风\\ &你是我遥不可及的梦\ ”\\ ——&\text{《你是我遥不可及的梦》 By 苍穹} \end{split} \]> Link 你是我遥不可及的梦-Bilibili
标签：Kingdom,frac,OI,卷积,sum,Assassination,Complex,const,pi 来源： https://www.cnblogs.com/LuckyGlass-blog/p/14379535.html