【Bluestein's Algorithm】[POJ2821]TN's Kingdom III - Assassination

题目大意

TN要暗杀Dzx，为了保密，他想到了这样一种方式：首先，把信息编码为N个实数，组成序列α，之后再随便搞一个长度为N的实数序列β。然后按照下面的步骤计算序列γ：
0、做一个空序列γ。
1、把β倒过来。
2、把β向右平移一个元素。最右侧的元素补到左边。
3、计算此时α和β对应元素的积的和。将其加到γ的末尾。
4、如果γ还不足N个元素，重复步骤2和3。

虽然这种加密方法是很弱的，可是那些笨蛋刺客们却没法破解。然后，这些东西就被Dzx拿到了，于是他就躲过了暗杀……现在作为历史学家的你得到了β和γ，你需要解出α来获得研究TN的材料。

分析

说得更直白一点，$\gamma$就是$\alpha$和$\beta$的循环卷积。

运用$FFT$所进行的卷积本身就是循环卷积

我们假设$c$为$a$、$b$两个序列的卷积，即$c_k=\sum_{i=0}^ka_k\times b_{k-i}$
令n为大于等于c的长度而且是2的整数次幂的数，那么卷积也可以写成这个形式$c_k=\sum_{i,j}a_kb_j[i+j \equiv 0\pmod n]$
FFT就是根据第二个式子算出来的，当$a,b$长度相等而且$n$等于他们的长度的时候，根据这个式子算出来的就是循环卷积。平时计算的时候由于$n\ge c$的长度，所以和线性卷积的式子没有区别。
很显然，我们对$\gamma$和$\beta$做$DFT$，然后相除，再$IDFT$即可，但是用$FFT$进行$DFT$要求长度必须是2的整数次幂，如果原长不是2的整数次幂，我们强行弄成2的整数次幂的话，在做$\gamma$的$DFT$和最后$IDFT$时就会出现问题，因为你根本不知道你这样做完了是一个什么东西。
所以，我们需要一个可以做任意长度卷积的东西，$Bluestein's\ \ Algorithm$或者混合基$FFT$。

Bluestein’s Algorithm

我们考虑$DFT$的式子

$$\begin{align}A_k&=\sum _{j=0}^{N-1}a_j\omega ^{jk}_N\\&=\sum _{j=0}^{N-1}a_je^{\frac{2\pi i}{N} jk} \\jk&=\frac{-(j-k)^2+j^2+k^2}{2}\\A_k&=\sum _{j=0}^{N-1}a_j e ^{\frac{\pi i}{N}(-(j-k)^2+j^2+k^2)}\\&=e^{\frac{\pi i}{N} k^2} \sum a_je^{\frac{\pi i}{N}j^2} e^{\frac{\pi i}{N}(-(k-j)^2)}\end{align}$$

设$X_j=a_je^{\frac{\pi i}{N}j^2},Y_j=e^{-\frac{\pi i}{N}j^2}$

$$A_k=e^{\frac{\pi i}{N}{k^2}}\sum _{j=0}^{N}X_jY_{k-j}$$
这是一个卷积的形式，我们发现我们可以通过一次卷积来做一次 $DFT$.
ps:其实是CZT，fft是dft的快速算法，其实就是N点dft算法，就是计算量小一点。N点dft的本质是z变换后，在z域单位圆上等间距N点连续采样。
$IDFT$的时候，可以类比用 $FFT$进行 $IDFT$时做出的改变即可。
我们发现 $k-j$可能小于0，我们只需要将 $Y(x)$右移 N N <script type="math/tex" id="MathJax-Element-58">N</script>位即可。

代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const double pi=acos(-1);
#define MAXN (1<<17)
int n;
struct cpx{
    double r,i;
    inline cpx(){
    }
    inline cpx(double r,double i=0):r(r),i(i){
    }
    inline cpx operator+(const cpx &b)const{
        return cpx(r+b.r,i+b.i);
    }
    inline cpx operator-(const cpx &b)const{
        return cpx(r-b.r,i-b.i);
    }
    inline cpx operator*(const cpx &b)const{
        return cpx(r*b.r-i*b.i,r*b.i+i*b.r);
    }
    inline cpx operator/(const cpx &b)const{
        return cpx((r*b.r+i*b.i)/(b.r*b.r+b.i*b.i),(i*b.r-r*b.i)/(b.r*b.r+b.i*b.i));
    }
    inline cpx& operator*=(const cpx &b){
        return *this=*this*b;
    }
    inline cpx &operator/=(const cpx &b){
        return *this=*this/b;
    }
}beta [MAXN+10],gamma[MAXN+10],alpha[MAXN+10],A[MAXN*4+10],B[MAXN*4+10];
void fft(cpx *a,int N,int f){
    int i,j,t,k;
    for(i=1,j=0;i<N-1;i++){
        for(t=N;j^=(t>>=1),~j&t;);
        if(i<j)
            swap(a[i],a[j]);
    }
    for(i=1;i<N;i<<=1){
        cpx wn(cos(pi/i),f*sin(pi/i));
        t=i<<1;
        for(j=0;j<N;j+=t){
            cpx w(1,0);
            for(k=0;k<i;k++,w*=wn){
                cpx x(a[j+k]),y(w*a[j+i+k]);
                a[j+k]=x+y;
                a[j+i+k]=x-y;
            }
        }
    }
    if(f==-1){
        for(i=0;i<N;i++)
            a[i]/=N;
    }
}
void bluestein(cpx *a,int n,int f){
    int N,i;
    memset(A,0,sizeof A);
    memset(B,0,sizeof B);
    for(i=0;i<n;i++)
        A[i]=cpx(cos(pi*i*i/n),f*sin(pi*i*i/n))*a[i];
    for(i=0;i<(n<<1);i++)
        B[i]=cpx(cos(pi*(i-n)*(i-n)/n),-f*sin(pi*(i-n)*(i-n)/n));
    for(N=1;N<(n<<2);N<<=1);
    fft(A,N,1);
    fft(B,N,1);
    for(i=0;i<N;i++)
        A[i]*=B[i];
    fft(A,N,-1);
    for(i=0;i<n;i++){
        a[i]=A[i+n]*cpx(cos(pi*i*i/n),f*sin(pi*i*i/n));
        if(f==-1)
            a[i]/=n;
    }
}
void read(){
    scanf("%d",&n);
    int i;
    for(i=0;i<n;i++)
        scanf("%lf",&beta[i].r);
    for(i=0;i<n;i++)
        scanf("%lf",&gamma[i].r);
}
void solve(){
    bluestein(beta,n,1);
    bluestein(gamma,n,1);
    for(int i=0;i<n;i++)
        alpha[i]=gamma[i]/beta[i];
    bluestein(alpha,n,-1);
}
void print(){
    int i;
    for(i=0;i<n;i++)
        printf("%.4f\n",alpha[i].r);
}
int main()
{
    read();
    solve();
    print();
}