BZOJ3527 力

题目

给出n个数qi，给出Fj的定义如下：

令Ei=Fi/qi，求Ei.

思路

年轻人的第一道FFT。

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

设：

\[f[i]=q_i\\ g[i]=\frac{1}{i^2} \]

那么原式可以化为：

\[E_j=\sum_{i<j}f[i]g[j-i]-\sum_{i>j}f[i]g[i-j] \]

也就是:

\[E_j=\sum_{i=0}^{j-1}f[i]g[j-i]-\sum_{i=j+1}^{n}f[i]g[i-j] \]

前半部分是卷积的形式，显然可以用FFT维护，重点是后面的部分。

我们来变一下形：

\[\sum_{i=j+1}^{n}f[i]g[i-j]\\ 令i=i-j\\ \sum_{i=0}^{n-j}f[i+j]g[i] \]

然后将\(i+j+i\)凑成\(n-j\),也就是令\(f'[n-j-i]=f[i+j]\)

再简单变形：\(f'[i]=f[n-i]\)

于是\(f'[]\)也构造出来了，这样两边都化为了卷积的形式，就可以用FFT求解了。

代码

#include<bits/stdc++.h>
#define M 400005
using namespace std;
int n;
struct Complex{
	double x,y;
	Complex(){}
	Complex(double _x,double _y):x(_x),y(_y){}
	Complex operator + (const Complex &res) const{
		return (Complex){x+res.x,y+res.y};	
	}
	Complex operator - (const Complex &res) const{
		return (Complex){x-res.x,y-res.y};	
	}
	Complex operator * (const Complex &res) const{
		return (Complex){x*res.x-y*res.y,x*res.y+y*res.x};	
	}
};
Complex A[M],B[M];
double pi=acos(-1.0);
void FFT(Complex *y,int n,int f){
	if(n==1)return;
	Complex L[n>>1],R[n>>1];
	for(int i=0;i<n;i+=2)L[i>>1]=y[i],R[i>>1]=y[i+1];
	FFT(L,n>>1,f);FFT(R,n>>1,f);
	Complex wn(cos(2*pi/n),f*sin(2*pi/n)),w(1,0);
	for(int i=0;i<(n>>1);i++,w=w*wn){
		y[i]=L[i]+w*R[i];
		y[i+(n>>1)]=L[i]-w*R[i];
	}
}
double q[M],b[M],C[M];
int main(){
	cin>>n;n--;
	for(int i=0;i<=n;i++)scanf("%lf",&q[i]);
	for(int i=0;i<=n;i++)A[i]=(Complex){q[i],0};
	for(int i=1;i<=n;i++){
		b[i]=1.0/i/i;
		B[i]=Complex(b[i],0);
	}
	int nn=n,mm=n;
	mm+=nn;
	for(nn=1;nn<=mm;nn<<=1);
	FFT(A,nn,1);FFT(B,nn,1);
	for(int i=0;i<=nn;i++)A[i]=A[i]*B[i];
	FFT(A,nn,-1);
	for(int i=0;i<=n;i++)C[i]=A[i].x/nn;
	for(int i=0;i<=nn;i++)A[i]=B[i]=Complex(0,0);
	for(int i=0;i<=n;i++)A[i]=Complex(q[n-i],0);
	for(int i=1;i<=n;i++)B[i]=Complex(b[i],0);
	FFT(A,nn,1);FFT(B,nn,1);
	for(int i=0;i<=nn;i++)A[i]=A[i]*B[i];
	FFT(A,nn,-1);
	for(int i=0;i<=n;i++)C[i]-=A[n-i].x/nn;
	for(int i=0;i<=n;i++)printf("%.3f\n",C[i]);
	return 0;
}