BZOJ3527: [Zjoi2014]力

【传送门：BZOJ3527】

---

简要题意：

　　给出n个数q[i]，给出F[j]的定义：

　　$$F[j]=\sum_{i<j}\frac{q[i]q[j]}{(i-j)^2}-\sum_{i>j}\frac{q[i]q[j]}{(i-j)^2}$$

　　令E[i]=F[i]/q[i]，求E[i]

---

题解：

　　最近刷FFT，心力憔悴

　　对于这道题直接求F[j]显然不理想

　　我们可以把原式转化为：

　　$$E[i]=\sum_{j<i}\frac{q[j]}{(i-j)^2}-\sum_{j>i}\frac{q[j]}{(i-j)^2}$$

　　为了使这个式子变为卷积的形式，我们发现$\sum_{j<i}$其实等效于$\sum_{j=0}^{i-1}$，而$\sum_{j>i}$其实也等效于$\sum_{j=i+1}^{n-1}$

　　所以我们将式子转化：

　　$$E[i]=\sum_{j=0}^{i-1}\frac{q[j]}{(i-j)^2}-\sum_{j=i+1}^{n-1}\frac{q[j]}{(i-j)^2}$$

　　设$f(i)=q[i]$,$g(i)=\frac{1}{i^2}$

　　所以式子可以化为：

　　$$E[i]=\sum_{j=0}^{i-1}f(j)g(i-j)-\sum_{j=i+1}^{n-1}f(j)g(i-j)$$

　　$$E[i]=\sum_{j=0}^{i-1}f(j)g(i-j)-\sum_{j=0}^{n-i-2}f(j-i-1)g(2*i-j+1)$$

　　然而这个式子求卷积有点小东西要处理，为了避免，我们把它变成这样：

　　$$E[i]=\sum_{j=0}^{i}f(j)g(i-j)-f(i)g(i-i)-\sum_{j=0}^{n-i-1}f(j+i)g(-j)+f(i)g(i-i)$$

　　$$E[i]=\sum_{j=0}^{i}f(j)g(i-j)-\sum_{j=0}^{n-i-1}f(j+i)g(j)$$

　　上面的右边的式子$g(j)$，本应是$g(-j)$，但是两者的值相等，所以还是统一用$g(j)$

　　左边的式子显然可以直接FFT求卷积得出，接下来只要处理右边的式子就可以了

　　对于式子$\sum_{j=0}^{n-i-1}f(j+i)g(j)$，显然不能直接求

　　所以我们设$f'(i)=f(n-i-1)$，那么上述式子就可以转化为:$\sum_{j=0}^{n-i-1}f'(n-i-j-1)g(j)$

　　设$X(i)=\sum_{j=0}^{n-i-1}f'(n-i-j-1)g(j)$，那么$X(n-i-1)=\sum_{j=0}^{i}f'(i-j)g(j)$

　　显然$X(n-i-1)=\sum_{j=0}^{i}f'(i-j)g(j)$可以直接用FFT求

　　那么答案就是：

　　$$E[i]=\sum_{j=0}^{i}f(j)g(i-j)-X(i)$$

---

参考代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
using namespace std;
const double PI=acos(-1.0);
struct Complex
{
    double r,i;
    Complex(){}
    Complex(double _r,double _i){r=_r;i=_i;}
    friend Complex operator + (const Complex &x,const Complex &y){return Complex(x.r+y.r,x.i+y.i);}
    friend Complex operator - (const Complex &x,const Complex &y){return Complex(x.r-y.r,x.i-y.i);}
    friend Complex operator * (const Complex &x,const Complex &y){return Complex(x.r*y.r-x.i*y.i,x.r*y.i+x.i*y.r);}
}a[410000],b[410000];
int n,m;
int R[410000];
void fft(Complex *y,int len,int on)
{
    for(int i=0;i<len;i++) if(i<R[i]) swap(y[i],y[R[i]]);
    for(int i=1;i<len;i<<=1)
    {
        Complex wn(cos(PI/i),sin(on*PI/i));
        for(int j=0;j<len;j+=(i<<1))
        {
            Complex w(1,0);
            for(int k=0;k<i;k++,w=w*wn)
            {
                Complex u=y[j+k];
                Complex v=w*y[j+k+i];
                y[j+k]=u+v;
                y[j+k+i]=u-v;
            }
        }
    }
    if(on==-1) for(int i=0;i<len;i++) y[i].r/=len;
}
void calc(int n)
{
    int L=0,m=2*n;
    for(n=1;n<=m;n<<=1) L++;
    memset(R,0,sizeof(R));
    for(int i=0;i<n;i++) R[i]=(R[i>>1]>>1)|(i&1)<<(L-1);
    fft(a,n,1);
    fft(b,n,1);
    for(int i=0;i<=n;i++) a[i]=a[i]*b[i];
    fft(a,n,-1);
}
double q[410000];
double d[410000];
int main()
{
    int n;
    scanf("%d",&n);n--;
    for(int i=0;i<=n;i++)
    {
        scanf("%lf",&q[i]);
        a[i].r=q[i];
    }
    for(int i=1;i<=n;i++) b[i].r=1.0/(double(i)*double(i));
    calc(n);
    memset(d,0,sizeof(d));
    for(int i=0;i<=n;i++) d[i]+=a[i].r;
    memset(a,0,sizeof(a));
    memset(b,0,sizeof(b));
    for(int i=0;i<=n;i++) a[i].r=q[n-i];
    for(int i=1;i<=n;i++) b[i].r=1.0/(double(i)*double(i));
    calc(n);
    for(int i=0;i<=n;i++) d[i]-=a[n-i].r;
    for(int i=0;i<=n;i++) printf("%.3lf\n",d[i]);
    return 0;
}