【jzoj3617】【ZJOI2014】【力】【fft】

来源：互联网发布：异形契约知乎编辑：程序博客网时间：2024/05/16 18:13

题目大意

解题思路

Ei=∑j<iqj/(i−j)2−∑j>iqj/(i−j)2

我们设f[i]=1/i/i，gi=∑j<iqj/(i−j)2=∑j<iqjfi−j

我们发现这个数一个卷积，我们把q和f当作两个从左到右是从低到高位的数做乘法（在最前面各加一个零位，消去当前位的影响），由于乘法的定义，乘后的第i位即为g[i]。对于后面的部分可以类似处理。

code

#include<cmath>#include<cstdio>#include<algorithm>#define LF double#define LL long long#define min(a,b) ((a<b)?a:b)#define max(a,b) ((a>b)?a:b)#define fo(i,j,k) for(int i=j;i<=k;i++)#define fd(i,j,k) for(int i=j;i>=k;i--)using namespace std;int const mxn=1e5;int n,N,M;LF pi=acos(-1);struct rec{    LF x,y;    rec(LF X=0,LF Y=0){x=X;y=Y;}};rec operator+(rec x,rec y){return rec(x.x+y.x,x.y+y.y);}rec operator-(rec x,rec y){return rec(x.x-y.x,x.y-y.y);}rec operator*(rec x,rec y){return rec(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);}rec a[mxn*4+10],b[mxn*4+10],c[mxn*4+10],d[mxn*4+10],t[mxn*4+10];void DFT(rec *a,int tag){    fo(i,0,N-1){        int pos=0;        for(int ii=i,bit=1;bit<=M;pos=(pos<<1)+(ii&1),ii=ii>>1,bit++);        t[pos]=a[i];    }    for(int size=2;size<=N;size=size<<1){        int half=size>>1;        fo(i,0,half-1){            rec w=rec(cos(tag*pi*i/half),sin(tag*pi*i/half));            for(int j=i;j<N;j+=size){                rec u=t[j],v=w*t[j+half];                t[j]=u+v;                t[j+half]=u-v;            }        }    }    fo(i,0,N-1)a[i]=t[i];}void FFT(rec *a,rec *b){    DFT(a,1);    fo(i,0,N-1)a[i]=a[i]*b[i];    DFT(a,-1);    fo(i,0,N-1)a[i].x=a[i].x/N;}int up(LF x){    return int(x)+(int(x)!=x);}int main(){    freopen("d.in","r",stdin);    freopen("d.out","w",stdout);    scanf("%d",&n);    M=up(log(n*2+2)/log(2)),N=pow(2,M);    fo(i,1,n)scanf("%lf",&a[i].x),b[n-i+1].x=a[i].x;    fo(i,1,n)c[i].x=1.0/i/i;    DFT(c,1);    FFT(a,c);    FFT(b,c);    fo(i,1,n)printf("%lf\n",a[i].x-b[n-i+1].x);    return 0;}

0 0