BZOJ 3527: [Zjoi2014]力

Description

求 \(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}\)

Sol

FFT.

我们可以发现他是一个卷积的形式,每次从\(i^2\) 卷到 \((n-i-1)^2\) .

既然是卷积,那么直接FFT就好了,但是FFT是让指数相等,也就是这里面的下标相等,所以必须要翻转这两个数组其中一个就可以了,随便翻就行.

然后从某一个下下标位置开始输出.

Code

/**************************************************************
    Problem: 3527
    User: BeiYu
    Language: C++
    Result: Accepted
    Time:2696 ms
    Memory:24732 kb
****************************************************************/
 
#include <bits/stdc++.h>
using namespace std;
 
#define mpr make_pair
#define rr first
#define ii second
const int N = 5e5+50;
const double Pi = M_PI;
typedef pair< double,double > Complex;
 
Complex operator + (const Complex &a,const Complex &b) {
    return mpr(a.rr+b.rr,a.ii+b.ii);
}
Complex operator - (const Complex &a,const Complex &b) {
    return mpr(a.rr-b.rr,a.ii-b.ii);
}
Complex operator * (const Complex &a,const Complex &b) {
    return mpr(a.rr*b.rr-a.ii*b.ii,a.rr*b.ii+a.ii*b.rr);
}
 
int n,m;
Complex a[N],b[N],c[N];
 
void init(int x) {
    for(n=1,m=0;n<x;n<<=1,m++);
    m++,n<<=1;
}
void Rev(Complex a[]) {
    for(int i=0,j=0,k;i<n;i++) {
        if(i>j) swap(a[i],a[j]);
        for(k=n>>1;(j^=k)<k;k>>=1);
    }
}
void DFT(Complex y[],int r) {
    Rev(y);
    for(int i=2;i<=n;i<<=1) {
        Complex wi=mpr(cos(2.0*Pi/i),r*sin(2.0*Pi/i));
        for(int k=0;k<n;k+=i) {
            Complex w=mpr(1.0,0.0);
            for(int j=k;j<k+i/2;j++) {
                Complex t1=y[j],t2=w*y[j+i/2];
                y[j]=t1+t2,y[j+i/2]=t1-t2;
                w=w*wi;
            }
        }
    }if(r==-1) for(int i=0;i<n;i++) y[i].rr/=n;
}
void FFT(Complex a[],Complex b[],Complex c[]) {
    DFT(a,1);DFT(b,1);
    for(int i=0;i<n;i++) c[i]=a[i]*b[i];
    DFT(c,-1);
}
int main() {
//  freopen("in.in","r",stdin);
    int l;
    scanf("%d",&l);
    for(int i=0;i<l;i++) scanf("%lf",&a[i].rr);
    for(int i=0;i<l-1;i++) {
        b[i]=mpr(1.0/(l-i-1.0)/(l-i-1.0),0.0);
        b[2*l-1-i-1]=mpr(-1.0/(l-i-1.0)/(l-i-1.0),0.0);
    }
    init(l);
//  for(int i=0;i<n;i++) cout<<b[i].rr<<" ";cout<<endl;
//  for(int i=0;i<n;i++) cout<<a[i].rr<<" ";cout<<endl;
//  reverse(a,a+l);
    reverse(b,b+2*l-1);
//  reverse(b,b+n);
//  for(int i=0;i<n;i++) cout<<b[i].rr<<" ";cout<<endl;
//  for(int i=0;i<n;i++) cout<<a[i].rr<<" ";cout<<endl;
    FFT(a,b,c);
//  for(int i=0;i<n;i++) printf("%lf\n",c[i].rr);
    for(int i=l-1;i<2*l-1;i++) printf("%lf\n",c[i].rr);
    return 0;
}
/*
5 
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880
*/