P3338 [ZJOI2014]力
蒟蒻第一次做不是板子(?不是吗)的fft题(竟然很流畅地50minA了)
首先化简式子
$E_{i}=\sum_{j<i}^{ }\frac{q_{j}}{(i-j)^{2}}-\sum_{j>i}^{ }\frac{q_{j}}{(i-j)^{2}}$
然后可以设数列{$q_{i}$}与{$a_{i}$}
{$a_{i}$}的定义可以看代码,不太好描述
然后$ans_{i}=\sum_{k+j=n+i}q_{j}*a_{k}$(不知道写法有没有问题...)
然后可以发现这是卷积的形式
然后fft就可以快乐卷积辣
code:
#include<bits/stdc++.h> using namespace std; const int N=1000000; const double pi=acos(-1); struct comp{ double x,y; comp(double xx=0,double yy=0){ x=xx,y=yy; } }; comp operator +(comp a,comp b){ return comp(a.x+b.x,a.y+b.y); } comp operator -(comp a,comp b){ return comp(a.x-b.x,a.y-b.y); } comp operator *(comp a,comp b){ return comp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x); } comp a[N],b[N]; int n,L,lim=1,r[N]; void fft(comp *f,int op){ for(int i=0;i<lim;i++){ if(i<r[i])swap(f[i],f[r[i]]); } for(int p=2;p<=lim;p<<=1){ int len=p>>1; comp tem(cos(pi/len),op*sin(pi/len)); for(int k=0;k<lim;k+=p){ comp buf(1,0); for(int j=k;j<k+len;j++){ comp tt=buf*f[j+len]; f[j+len]=f[j]-tt; f[j]=f[j]+tt; buf=buf*tem; } } } } int main(){ cin>>n; for(int i=0;i<n;i++){ double x; scanf("%lf",&x); b[i].x=x; } for(int i=-n+1;i<=n-1;i++){ if(i<0) a[i+n-1].x=-1.0/i/i; else if(i==0) a[i+n-1].x=0; else a[i+n-1].x=1.0/i/i; } while((1<<L)<(3*n-3))L++,lim<<=1; for(int i=0;i<lim;i++){ r[i]=(r[i>>1]>>1)|((i&1)<<(L-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("%.3f\n",a[i].x/lim); } return 0; }