LG P3338 [ZJOI2014]力

Description

给出 $n$ 个数 $q_1,q_2, dots q_n$，定义

$$F_j~=~sum_{i = 1}^{j - 1} frac{q_i imes q_j}{(i - j)^2}~-~sum_{i = j + 1}^{n} frac{q_i imes q_j}{(i - j)^2}$$

$$E_i~=~frac{F_i}{q_i}$$

对 $1 leq i leq n$，求 $E_i$ 的值。

Solution

原式可化为$$sum_{i=1}^{j-1} frac{q_i}{(i-j)^2} - sum_{i=j+1}^n frac{q_i}{(i-j)^2}$$

设$f_i=q_i,g_i=frac 1{i^2}$

令$f_0=g_0=0$

则原式化为$$sum_{i=1}^j f_i imes g_{j-i} - sum_j^n f_i imes g_{i-j}$$

设$f'_i = f_{n-i}$，则原式化为

$$sum_{i=0}^j f_i imes g_{j-i} - sum_{i=0}^{n-j} f'_{n-j-i} imes g_i$$

两次FFT求解

#include<iostream>
#include<complex>
#include<cstdio>
#include<cmath>
using namespace std;
int n,s=2,tot=1,rev[400005];
const double pi=acos(-1);
complex<double> f[400005],g[400005],ff[400005];
inline int read()
{
    int f=1,w=0;
    char ch=0;
    while(ch<'0'||ch>'9'){if(ch=='-') f=-1;ch=getchar();}
    while(ch>='0'&&ch<='9') w=(w<<1)+(w<<3)+ch-'0',ch=getchar();
    return f*w;
}
void fft(complex<double>*a,int n,int inv)
{
    for(int i=0;i<n;i++) if(i<rev[i]) swap(a[i],a[rev[i]]);
    for(int i=1;i<n;i*=2)
    {
        complex<double>wn=exp(complex<double>(0,inv*pi/i));
        for(int j=0;j<n;j+=i*2)
        {
            complex<double>w(1,0);
            for(int k=j;k<i+j;k++)
            {
                complex<double>x=a[k],y=w*a[k+i];
                a[k]=x+y,a[k+i]=x-y,w*=wn;
            }
        }
    }
    if(inv==-1) for(int i=0;i<n;i++) a[i]/=n;
}
int main()
{
    n=read();
    for(int i=1;i<=n;i++) scanf("%lf",&f[i]),ff[n-i]=f[i];
    for(int i=1;i<=n;i++) g[i]=1.0/i/i;
    while(s<n*2) s<<=1,tot++;
    for(int i=0;i<s;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(tot-1));
    fft(f,s,1),fft(g,s,1),fft(ff,s,1);
    for(int i=0;i<s;i++) f[i]*=g[i];
    for(int i=0;i<s;i++) ff[i]*=g[i];
    fft(f,s,-1),fft(ff,s,-1);
    for(int i=1;i<=n;i++) printf("%.3lf
",f[i].real()-ff[n-i].real());
    return 0;
}