BZOJ 3527 ZJOI2014 力 快速傅里叶变换

题目大意：给定n个点，第i个点和第j个点之间的库仑力为(qi*qj)/(i-j)^2，定义左侧为正方向，求每个点受的合力与电荷量的比值

题解详见 http://eolv.farbox.com/post/shui-yu-zheng-feng/2014-12-07 我懒得打了- -

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 263000
#define PI 3.1415926535897932384626433832795028841971
using namespace std;
struct Complex{
	long double a,b;
	Complex() {}
	Complex(long double _,long double __):a(_),b(__) {}
	Complex operator + (const Complex &x) const
	{
		return Complex(a+x.a,b+x.b);
	}
	Complex operator - (const Complex &x) const
	{
		return Complex(a-x.a,b-x.b);
	}
	Complex operator * (const Complex &x) const
	{
		return Complex(a*x.a-b*x.b,a*x.b+b*x.a);
	}
	void operator *= (const Complex &x)
	{
		*this=(*this)*x;
	}
}a[M],b[M],c[M];
int n;
long double q[M],ans[M];
void FFT(Complex x[],int n,int type)
{
	static Complex temp[M];
	if(n==1) return ;
	int i;
	for(i=0;i<n;i+=2)
		temp[i>>1]=x[i],temp[i+n>>1]=x[i+1];
	memcpy(x,temp,sizeof(Complex)*n);
	Complex *l=x,*r=x+(n>>1);
	FFT(l,n>>1,type);FFT(r,n>>1,type);
	Complex root(cos(type*2*PI/n),sin(type*2*PI/n)),w(1,0);
	for(i=0;i<n>>1;i++,w*=root)
		temp[i]=l[i]+w*r[i],temp[i+(n>>1)]=l[i]-w*r[i];
	memcpy(x,temp,sizeof(Complex)*n);
}
int main()
{
	int i,digit;
	double x;
	
	cin>>n;
	for(digit=1;digit<=n<<1;digit<<=1);

	for(i=0;i<n;i++)
		scanf("%lf",&x),q[i]=x;

	for(i=1;i<n;i++)
		b[i].a=1.0/i/i;
	FFT(b,digit,1);
	
	for(i=0;i<n;i++)
		a[i]=Complex(q[i],0);
	FFT(a,digit,1);
	for(i=0;i<digit;i++)
		c[i]=a[i]*b[i];
	FFT(c,digit,-1);
	for(i=0;i<n;i++)
		ans[i]+=c[i].a;

	for(i=0;i<n;i++)
		a[i]=Complex(q[n-i-1],0);
	for(i=n;i<digit;i++)
		a[i]=Complex(0,0);
	FFT(a,digit,1);
	for(i=0;i<digit;i++)
		c[i]=a[i]*b[i];
	FFT(c,digit,-1);
	for(i=0;i<n;i++)
		ans[i]-=c[n-i-1].a;

	for(i=0;i<n;i++)
		printf("%lf\n",(double)ans[i]/digit);
		//cout<<ans[i]/digit<<endl;
	return 0;
}