#5 #9 77分,迭代FFT求助
查看原帖
#5 #9 77分,迭代FFT求助
616090
hahaxy3楼主2024/8/10 19:01
#include <bits/stdc++.h>
const double pi=acos(-1);
using namespace std;
struct acomplex{
    double real=0.0,imag=0.0;
};
vector<acomplex> a3,b3;
acomplex operator + (acomplex a1,acomplex b1){
    acomplex c;
    c.real=a1.real+b1.real;
    c.imag=a1.imag+b1.imag;
    return c;
}
acomplex operator - (acomplex a1,acomplex b1){
    acomplex c;
    c.real=a1.real-b1.real;
    c.imag=a1.imag-b1.imag;
    return c;
}
acomplex operator * (acomplex a1,acomplex b1){
    acomplex c;
    c.real=a1.real*b1.real-a1.imag*b1.imag;
    c.imag=a1.real*b1.imag+b1.real*a1.imag;
    return c;
}
int lim=1,n=0,o,p;
vector<int>ind,ak;
//vector<acomplex>a3,b3;
void fft(vector<acomplex> &a4,int opt){
    for(int i=0;i<lim;i++){
    if(ind[i]>i)swap(a4[i],a4[ind[i]]);
    }
    for(int i=1;i<=n;i++){
    int n1=1<<i;
    acomplex wn;
    wn.real=cos(2.0*pi/n1);
    wn.imag=sin(2.0*pi/n1)*opt;
    for(int j=0;j<lim;j+=n1){
        acomplex w;
        w.real=1.0;
        w.imag=0.0;
        for(int k=0;k<n1/2;k++){
        //cout<<j+k<<" ";
        acomplex u=a4[j+k];
        acomplex v=w*a4[j+k+n1/2];
        a4[j+k]=u+v;
        a4[j+k+n1/2]=u-v;
        w=w*wn;
        }
    }
    }
    if(opt==-1)for(int i=0;i<lim;i++)a4[i].real/=lim;
}
int main(){
    ios::sync_with_stdio(0);
    cin>>o>>p;
    o++;
    p++;
    while(lim<(o+p)){
    lim<<=1;
    n++;
    }
    a3.resize(lim);
    b3.resize(lim);
    for(int i=0;i<o;i++)cin>>a3[i].real;
    for(int i=0;i<p;i++)cin>>b3[i].real;
    ind.push_back(0);
    ind.push_back(1<<(n-1));
    for(int inde=1;inde<lim/2;inde++){
    ind.push_back(ind[inde]>>1);
    ind.push_back((ind[inde]>>1)|(1<<(n-1)));
    }
    //for(int i=0;i<lim;i++)cout<<ind[i]<<" ";
    //a3.resize(lim),b3.resize(lim);
    fft(a3,1);fft(b3,1);
    for(int i=0;i<lim;i++){
    a3[i]=a3[i]*b3[i];
    //cout<<a3[i].real<<'+'<<a3[i].imag<<'i'<<" ";
    }
    fft(a3,-1);
    ak.resize(lim);
    for(int i=0;i<lim;i++){
    ak[i]=((int)(a3[i].real+0.5));
    //ak[i+1]+=(((int)(a3[i]).real+0.5))/10;
    }
    for(int i=0;i<lim;i++)if(ak[i])cout<<ak[i]<<" ";
    return 0;
}
2024/8/10 19:01
加载中...