[bzoj2194] fast Fourier II FFT

Keywords: less

Main idea:

Please calculate C[k]=sigma(a[i]*b[i-k]), where k < = I < n, and N < = 10 ^ 5. a. The elements in B are non negative integers less than or equal to 100.

Train of thought:

It is not difficult to find that the subscript difference of a pair of contributing points in two arrays is a fixed value. A classic way is to reverse one of the arrays, so that the sum becomes a fixed value, which just satisfies the coefficient representation of polynomial multiplication. Direct FFT is enough.

#include<bits/stdc++.h>

#define REP(i,a,b) for(int i=a,i##_end_=b;i<=i##_end_;++i)
#define DREP(i,a,b) for(int i=a,i##_end_=b;i>=i##_end_;--i)
typedef long long ll;

using namespace std;

void File(){
    freopen("bzoj2194.in","r",stdin);
    freopen("bzoj2194.out","w",stdout);
}

const double pi=acos(-1.0);
const int maxn=4e5+10;
int n,lim=1,len,dn[maxn],ans[maxn];

struct Complex{
    double x,y;
    Complex(double xx=0,double yy=0){x=xx,y=yy;}
};
Complex operator + (Complex a,Complex b){return (Complex){a.x+b.x,a.y+b.y};}
Complex operator - (Complex a,Complex b){return (Complex){a.x-b.x,a.y-b.y};}
Complex operator * (Complex a,Complex b){return (Complex){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
Complex a[maxn],b[maxn];

void fft(Complex *A,int ty){
    REP(i,0,lim-1)if(i<dn[i])swap(A[i],A[dn[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        Complex w(cos(2.0*pi/(mid<<1)),ty*sin(2.0*pi/(mid<<1)));
        for(int L=0;L<lim;L+=mid<<1){
            Complex wk(1,0);
            REP(i,L,L+mid-1){
                Complex u=A[i],v=wk*A[i+mid];
                A[i]=u+v;
                A[i+mid]=u-v;
                wk=wk*w;
            }
        }
    }
}

int main(){
    File();
    scanf("%d",&n);
    int u,v;
    REP(i,1,n){
        scanf("%d%d",&u,&v);
        a[i-1].x=u;
        b[n-i].x=v;
    }
    while(lim<=((n-1)<<1))lim<<=1,++len;
    if(!len)len=1;
    REP(i,0,lim-1)dn[i]=(dn[i>>1]>>1)|((i&1)<<(len-1));
    fft(a,1);
    fft(b,1);
    REP(i,0,lim-1)a[i]=a[i]*b[i];
    fft(a,-1);
    REP(i,0,n-1)ans[i]=(int)(a[n+i-1].x/lim+0.5);
    REP(i,0,n-1)printf("%d\n",ans[i]);
    return 0;
}

Posted by glansing on Thu, 02 Jan 2020 17:30:54 -0800