ZJOI2014 力

一道简单的\(FFT\)题
题目链接

题意简述

给定一个公式\(E_i=\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}\)
求\(E\).

---

解析

先把公式抄下来\(E_i=\sum_{j<i}\frac{q_j}{(i-j)^2}-\sum_{j>i}\frac{q_j}{(i-j)^2}\)
我们令
\(A_i=q_i\)
\(i<0,B_i=-\frac{1}{i^2}\)
\(i=0,B_i=0\)
\(i>0,B_i=\frac{1}{i^2}\)
于是发现\(E_i=\sum_{j=0}^{n-1}A_j{B_{i-j}}\)
然后就是\(FFT\)板子题了.
注意读入的时候要平移一下,把\(B_i\)的\(i\)变成自然数.

代码如下

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<vector>
#define N (600010)
#define inf (0x7f7f7f7f)
#define rg register int
#define Label puts("NAIVE")
#define spa print(' ')
#define ent print('\n')
#define rand() (((rand())<<(15))^(rand()))
typedef long double ld;
typedef long long LL;
typedef unsigned long long ull;
using namespace std;
const ld PI=3.14159265359;
struct com{
	ld a,b;
	com(){a=b=0;}
	com operator +(com x){
		com res;
		res.a=a+x.a,res.b=b+x.b;
		return res;
	}
	com operator -(com x){
		com res;
		res.a=a-x.a,res.b=b-x.b;
		return res;
	}
	com operator *(com x){
		com res;
		res.a=a*x.a-b*x.b,res.b=a*x.b+b*x.a;
		return res;
	}
}a[N],b[N];
int n,rev[N],Lim,len;
void FFT(com *a,int tp){
	for(int i=0;i<Lim;i++)
	if(i<rev[i])swap(a[i],a[rev[i]]);
	for(int i=1;i<Lim;i<<=1){
		com w; w.a=cos(PI/i),w.b=tp*sin(PI/i);
		for(int R=i<<1,j=0;j<Lim;j+=R){
			com p; p.a=1;
			for(int k=j;k<j+i;k++,p=p*w){
				com x=a[k],y=p*a[k+i];
				a[k]=x+y,a[k+i]=x-y;
			}
		}
	}
	if(tp==-1)
	for(int i=0;i<Lim;i++)a[i].a=a[i].a/(ld)Lim;
}
int main(){
	scanf("%d",&n);
	for(int i=0;i<n;i++)
	scanf("%Lf",&a[i].a);
	for(int i=-n+1;i<=n-1;i++){
		if(i==0)b[i+n-1].a=0;
		else b[i+n-1].a=1.0/(ld)i/(ld)i*((i<0)?(-1):1);
	}
	for(Lim=1;Lim<=(n*3);Lim<<=1)len++;
	for(int i=0;i<Lim;i++)
	rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	FFT(a,1),FFT(b,1);
	for(int i=0;i<Lim;i++)
	a[i]=a[i]*b[i];
	FFT(a,-1);
	for(int i=n-1;i<=2*n-2;i++)
	printf("%.10Lf\n",a[i].a);
}