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('
')
#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
",a[i].a);
}