[Zjoi2014]力(FFT/NTT)

Posted 2021-03-20 chasedeath

[Zjoi2014]力(FFT，卷积)

题意：给定(n)个点电荷，排在单位数轴上，求每个点的场强

考虑每个(i)对于每个(j)的贡献，分析式子

(E=cfrac{q_i}{(j-i)^2})

令(f(x)=sum q_ix^i)

(g(x)=sum a_ix^i,a_i=i<0?-frac{1}{i^2}:frac{1}{i^2})

(g(x))每一项(x)的指数其实是(j-i)的值

求(f(x)cdot g(x))即可，注意负数系数的话偏移一下

是一个简单的作差卷积，适合作为(FFT)入门题

#include<bits/stdc++.h>
using namespace std;

#define double long double

#define reg register
typedef long long ll;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); } 
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); } 
#define double long double

char IO;
int rd(){
    int s=0,f=0;
    while(!isdigit(IO=getchar())) if(IO=='-') f=1;
    do s=(s<<1)+(s<<3)+(IO^'0');
    while(isdigit(IO=getchar()));
    return f?-s:s;
}

const double PI=acos(-1);

const int N=(1<<19)+4,P=998244353;
const int g=3;

bool be;
int n,m;
struct Cp{
    double x,y;
    Cp(){}
    Cp(double _x,double _y){ x=_x,y=_y; }
    Cp operator + (const Cp t){ return Cp(x+t.x,y+t.y); }
    Cp operator - (const Cp t){ return Cp(x-t.x,y-t.y); }
    Cp operator * (const Cp t){ return Cp(x*t.x-y*t.y,x*t.y+y*t.x); }
} a[N],b[N];
int rev[N];

void FFT(int n,Cp *a,int f) {
    rep(i,0,n-1) if(rev[i]>i) swap(a[i],a[rev[i]]);
    for(reg int i=1;i<n;i<<=1) {
        Cp w(cos(PI/i),f*sin(PI/i));
        for(reg int l=0;l<n;l+=i*2) {
            Cp e(1,0);
            for(reg int j=l;j<l+i;++j,e=e*w) {
                Cp t=a[j+i]*e;
                a[j+i]=a[j]-t;
                a[j]=a[j]+t;
            }
        }
    }
    if(f==-1) rep(i,0,n-1) a[i].x/=n;
}


bool ed;
int main(){
    n=rd();
    int R=1,c=-1;
    while(R<=n*3) R<<=1,c++;
    rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<c);
    rep(i,0,n-1) scanf("%Lf",&a[i].x);
    rep(i,-n+1,n-1) if(i) {
        b[i+n].x=1.0/(1.0*i*i);
        if(i<0) b[i+n].x*=-1;
    }
    FFT(R,a,1),FFT(R,b,1);
    rep(i,0,R) a[i]=a[i]*b[i];
    FFT(R,a,-1);
    rep(i,0,n-1) printf("%.3Lf
",a[i+n].x);
}