BZOJ3527: [Zjoi2014]力【FFT】

3527: [Zjoi2014]力

先除去qiq_iqi​，令g(i)=1i2g(i)=\frac{1}{i^2}g(i)=i21​式子变成Ei=∑j=1i−1qi∗g(i−j)−∑j=i+1nqi∗g(j−i)E_i=\sum_{j=1}^{i-1} q_i*g(i-j)-\sum_{j=i+1}^nq_i*g(j-i)Ei​=∑j=1i−1​qi​∗g(i−j)−∑j=i+1n​qi​∗g(j−i)

然后将后一半翻转Ei=∑j=1i−1qi∗g(i−j)−∑j=1n−iqn−i∗g(j−i)E_i=\sum_{j=1}^{i-1} q_i*g(i-j)-\sum_{j=1}^{n-i}q_{n-i}*g(j-i)Ei​=∑j=1i−1​qi​∗g(i−j)−∑j=1n−i​qn−i​∗g(j−i)

就是两个卷积相减辣，直接FFT

#include<cmath>
#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN=200005;
const double PI=acos(-1.0);
int n,Len,lg2,rev[MAXN<<1];
struct CP{
	double x,y;
	CP(double X=0,double Y=0){x=X,y=Y;}
	CP operator +(const CP b)const{return CP(x+b.x,y+b.y);}
	CP operator -(const CP b)const{return CP(x-b.x,y-b.y);}
	CP operator *(const CP b)const{return CP(x*b.x-y*b.y,x*b.y+y*b.x);} 
};
CP G[MAXN<<1],B1[MAXN<<1],B2[MAXN<<1],Ans1[MAXN<<1],Ans2[MAXN<<1];
void FFT(CP *A,int opt){
	for(int i=0;i<Len;i++) if(rev[i]<i) swap(A[i],A[rev[i]]);
	for(int i=1;i<Len;i<<=1){
		CP WN(cos(PI/i),opt*sin(PI/i));
		for(int j=0;j<Len;j+=i<<1){
			CP WNK(1,0);
			for(int k=0;k<i;k++,WNK=WN*WNK){
				CP x=A[j+k],y=WNK*A[j+k+i];
				A[j+k]=x+y;A[j+k+i]=x-y; 
			}
		}
	}
	if(opt==-1) for(int i=0;i<Len;i++) A[i].x/=Len;
}
void Getrev(){for(int i=0;i<Len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg2-1));}
int main(){
	scanf("%d",&n);n--;
	for(int i=0;i<=n;i++) scanf("%lf",&B1[i].x),B2[n-i].x=B1[i].x;
	for(int i=1;i<=n;i++) G[i].x=1.0/i/i;
	for(Len=1,lg2=0;Len<=n<<1;Len<<=1,lg2++);Getrev();
	FFT(G,1),FFT(B1,1),FFT(B2,1);
	for(int i=0;i<Len;i++) Ans1[i]=G[i]*B1[i];
	for(int i=0;i<Len;i++) Ans2[i]=G[i]*B2[i];
	FFT(Ans1,-1),FFT(Ans2,-1);
	for(int i=0;i<=n;i++) printf("%.8lf\n",Ans1[i].x-Ans2[n-i].x);
	return 0;
} 

转载于:https://www.cnblogs.com/XSamsara/p/10547942.html