P3338 [ZJOI2014]力

Posted 2021-02-13 dreagonm

思路

颓柿子的题目

要求求这样的一个式子
[ F_j=sum_{i<j}frac{q_iq_j}{(i-j)^2}-sum_{i>j}frac{q_iq_j}{(i-j)^2} ]
令(E_i=frac{F_i}{q_i})，求所有的(E_i)

对于Ei，显然可以
[ E_i=sum_{j=0}^{i-1}frac{q_j}{(i-j)^2}-sum_{j=i+1}^nfrac{q_j}{(i-j)^2} ]
前后没什么关联，可以分开考虑，首先考虑前面部分
[ sum_{j=0}^{i-1}frac{q_j}{(i-j)^2} ]
设(f(i)=q_i)，(g(i)=frac{1}{i^2})

则前面一部分变为
[ sum_{j=0}^{i-1}f(j)g(i-j) ]
变成了卷积的形式

后面一部分是
[ sum_{j=i+1}^{n}frac{q_j}{(i-j)^2} ]
同理变成
[ sum_{j=i+1}^nf(j)g(j-i) ]
似乎没什么办法，但是我们可以把(f)整个反过来（就是(f(1))与(f(n))交换，与(f(2))与(f(n-1))交换）

设交换之后的(f)为(f')，(f'_i=f_{n-i+1})

则有
[ sum_{j=0}^{i-1}f'(j)g(j-i) ]
则原来的第i项和反转后的第(n-i+1)项相同

对应相减即可

代码

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <cmath>
#define int long long
using namespace std;
const double pi = acos(-1.0);
struct complex{//a+bi
    double a,b;
    complex operator + (complex &bx){
        return (complex){a+bx.a,b+bx.b};
    }
    complex operator - (complex &bx){
        return (complex){a-bx.a,b-bx.b};
    }
    complex operator * (complex &bx){
        return (complex){a*bx.a-b*bx.b,b*bx.a+a*bx.b};
    }
    complex conj(void){
        return (complex){a,-b};
    }
};
int n;
double q[401000],ans[401000];
complex inv[401000],wnk[401000],a[401000],b[401000];
void init(int len){
    for(int i=0;i<len;i++){
        wnk[i]=(complex){cos(2*pi*i/len),sin(2*pi*i/len)};
        inv[i]=wnk[i].conj();
    }
}
void FFT(complex *a,complex *opt,int n){
    int lim=0;
    while((1<<lim)<n){
        lim++;
    }
    n=(1<<lim);
    for(int i=0;i<n;i++){
        int t=0;
        for(int j=0;j<lim;j++)
            if((i>>j)&1)
                t|=(1<<(lim-j-1));
        if(t<i)
            swap(a[i],a[t]);
    }
    for(int i=2;i<=n;i<<=1){
        int len=i/2;
        for(complex *j=a;j!=a+n;j+=i){
            for(int k=0;k<len;k++){
                complex t=j[k+len]*opt[n/i*k];
                j[k+len]=j[k]-t;
                j[k]=j[k]+t;
            }
        }
    }
}
void Do_FFT(int n){
    int lim=0;
    while((1<<lim)<n){
        lim++;
    }
    n=(1<<lim);
    init(n);
    FFT(a,wnk,n);
    FFT(b,wnk,n);
    for(int i=0;i<n;i++)
        a[i]=a[i]*b[i];
    FFT(a,inv,n);
    for(int i=0;i<n;i++)
        a[i].a/=n;
}
signed main(){
    // freopen("1.in","r",stdin);
    // freopen("test.out","w",stdout);
    scanf("%lld",&n);
    for(int i=1;i<=n;i++){
        scanf("%lf",&q[i]);
    }
    for(int i=1;i<=n;i++){
        a[i].a=q[i];
        b[i].a=1.0/(i*i);
        a[i].b=0;
        b[i].b=0;
    }
    Do_FFT(2*n+4);
    for(int i=1;i<=n;i++)
        ans[i]=a[i].a;
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    for(int i=1;i<=n;i++){
        a[i].a=q[i];
        b[i].a=1.0/(i*i);
        a[i].b=0;
        b[i].b=0;
    }
    for(int i=1,j=n;i<j;i++,j--)
        swap(a[i],a[j]);
    Do_FFT(2*n+4);
    for(int i=1;i<=n;i++)
        ans[i]-=a[n-i+1].a;
    for(int i=1;i<=n;i++)
        printf("%.3lf
",ans[i]);
    return 0;
}