P3338 [ZJOI2014]力

思路

颓柿子的题目

要求求这样的一个式子
\[
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}^n\frac{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\n",ans[i]);
    return 0;
}

原文地址：https://www.cnblogs.com/dreagonm/p/10559235.html