题面

给出两个长度为\(n\)的序列\(B,C\),已知\(A\)和\(B\)的循环卷积为\(C\),求\(A\).

\(n<2^{17}\)

分析

Bluestein算法的模板题，可以参考这篇博客

再探快速傅里叶变换(FFT)学习笔记(其三)(循环卷积的Bluestein算法+分治FFT+FFT的优化+任意模数NTT)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define maxn (1<<17)
const double pi=acos(-1.0);
using namespace std;
struct com{
    double real;
    double imag;
    com(){

    }
    com(double _real,double _imag){
        real=_real;
        imag=_imag;
    }
    com(double x){
        real=x;
        imag=0;
    }
    void operator = (const com x){
        this->real=x.real;
        this->imag=x.imag;
    }
    void operator = (const double x){
        this->real=x;
        this->imag=0;
    }
    friend com operator + (com p,com q){
        return com(p.real+q.real,p.imag+q.imag);
    }
    friend com operator + (com p,double q){
        return com(p.real+q,p.imag);
    }
    void operator += (com q){
        *this=*this+q;
    }
    void operator += (double q){
        *this=*this+q;
    }
    friend com operator - (com p,com q){
        return com(p.real-q.real,p.imag-q.imag);
    }
    friend com operator - (com p,double q){
        return com(p.real-q,p.imag);
    }
    void operator -= (com q){
        *this=*this-q;
    }
    void operator -= (double q){
        *this=*this-q;
    }
    friend com operator * (com p,com q){
        return com(p.real*q.real-p.imag*q.imag,p.real*q.imag+p.imag*q.real);
    }
    friend com operator * (com p,double q){
        return com(p.real*q,p.imag*q);
    }
    void operator *= (com q){
        *this=(*this)*q;
    }
    void operator *= (double q){
        *this=(*this)*q;
    }
    friend com operator / (com p,double q){
        return com(p.real/q,p.imag/q);
    }
    void operator /= (double q){
        *this=(*this)/q;
    }
    friend com operator / (com p,com q){//复数的除法,类似解二元一次方程,代入复数乘法公式解出答案
        return com((p.real*q.real+p.imag*q.imag)/(q.real*q.real+q.imag*q.imag),(p.imag*q.real-p.real*q.imag)/(q.real*q.real+q.imag*q.imag));
    }
    void print(){
        printf("%lf + %lf i ",real,imag);
    }
};

void fft(com *x,int *rev,int n,int type){
    for(int i=0;i<n;i++) if(i<rev[i]) swap(x[i],x[rev[i]]);
    for(int len=1;len<n;len*=2){
        int sz=len*2;
        com wn1=com(cos(2*pi/sz),type*sin(2*pi/sz));
        for(int l=0;l<n;l+=sz){
            int r=l+len-1;
            com wnk=1;
            for(int i=l;i<=r;i++){
                com tmp=x[i+len];
                x[i+len]=x[i]-wnk*tmp;
                x[i]=x[i]+wnk*tmp;
                wnk=wnk*wn1;
            }
        }
    }
    if(type==-1) for(int i=0;i<n;i++) x[i]/=n;
}
void bluestein(com *a,int n,int type){
    static com x[maxn*4+5],y[maxn*4+5];
    static int rev[maxn*4+5];
    memset(x,0,sizeof(x));
    memset(y,0,sizeof(y));
    int N=1,L=0;
    while(N<n*4){
        L++;
        N*=2;
    }
    for(int i=0;i<N;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
    for(int i=0;i<n;i++) x[i]=com(cos(pi*i*i/n),type*sin(pi*i*i/n))*a[i];
    for(int i=0;i<n*2;i++) y[i]=com(cos(pi*(i-n)*(i-n)/n),-type*sin(pi*(i-n)*(i-n)/n));
    fft(x,rev,N,1);
    fft(y,rev,N,1);
    for(int i=0;i<N;i++) x[i]*=y[i];
    fft(x,rev,N,-1);
    for(int i=0;i<n;i++){
        a[i]=x[i+n]*com(cos(pi*i*i/n),type*sin(pi*i*i/n));
        if(type==-1) a[i]/=n;//一定记得除以n,因为做一次Bluestein相当于一次FFT,IFFT最后要除n,这里也要除n
    }
}
void div(com *a,com *b,com *c,int n){//求解A*B=C
    bluestein(b,n,1);
    bluestein(c,n,1);
    for(int i=0;i<n;i++) a[i]=c[i]/b[i];
    bluestein(a,n,-1);
}

int n;
com a[maxn+5],b[maxn+5],c[maxn+5];
int main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%lf",&b[i].real);
    for(int i=0;i<n;i++) scanf("%lf",&c[i].real);
    div(a,b,c,n);
    for(int i=0;i<n;i++) printf("%.4f\n",a[i].real);
}