[zjoi2014] 力

传送门

Description

给出 (n)个数 (q_i)，给出 (E_j)的定义如下：

[E_j=sum_{i<j}frac{q_i}{(i-j)^2}-sum_{i>j}frac{q_i}{(i-j)^2} ]

求(E_i)

Solution 

前后两个式子分别化成卷积的形式

用(fft)来算

Code 

#include<bits/stdc++.h>
using namespace std;
#define db double
#define reg register

const int MN=3e5+5;

struct cp{
    db x,y;
    cp(db x=0,db y=0):x(x),y(y){}
    cp operator+(cp o)const{return cp(x+o.x,y+o.y);}
    cp operator-(cp o)const{return cp(x-o.x,y-o.y);}
    cp operator*(cp o)const{return cp(x*o.x-y*o.y,x*o.y+y*o.x);}
    inline void swap(cp& o){cp t=*this;*this=o;o=t;} //wrong af
}A[MN],B[MN],C[MN];

const db Pi=acos(-1.);

int pos[MN];

void fft(cp *a,int ty,int N)
{
    reg int i,j,k;cp w,wn,x,y;
    for(i=0;i<N;++i)if(i<pos[i])a[i].swap(a[pos[i]]);
    for(i=1;i<N;i<<=1)
    {
        wn=cp(cos(Pi/i),sin(Pi/i)*ty);
        for(j=0;j<N;j+=(i<<1))
            for(w=cp(1,0),k=0;k<i;++k,w=w*wn)
            {
                x=a[j+k];y=a[j+i+k]*w;
                a[j+k]=x+y;a[j+i+k]=x-y;
            }
    }
    if(ty==-1) for(i=0;i<N;++i) a[i].x/=N;
}

int main()
{
    reg int n,m,i;
    scanf("%d",&n);
    for(m=1;m<n;m<<=1);m<<=1;
    for(i=0;i<m;++i) pos[i]=(pos[i>>1]>>1)|((i&1)*(m>>1));
    
    for(i=1;i<=n;++i) C[i].x=1./i/i;
    for(i=1;i<=n;++i) scanf("%lf",&A[i].x),B[n-i+1].x=A[i].x;

    fft(A,1,m);fft(B,1,m);fft(C,1,m);
    for(i=0;i<m;++i) A[i]=A[i]*C[i],B[i]=B[i]*C[i];
    fft(A,-1,m);fft(B,-1,m);
    
    for(i=1;i<=n;++i) printf("%.10lf
",A[i].x-B[n-i+1].x);
    return 0;
}

---

Blog来自PaperCloud，未经允许，请勿转载，TKS！