題目描述

輸入輸出格式

輸入格式：

第一行一個整數n。

接下來n行每行輸入一個數，第i行表示qi。

輸出格式:

n行，第i行輸出Ei。

與標準答案誤差不超過1e-2即可。

輸入輸出樣例

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

說明

對於30%的數據，n≤1000。

對於50%的數據，n≤60000。

對於100%的數據，n≤100000，0<qi<1000000000。

題解：

E[i]=∑q[j]/(i-j)^2-∑q[j]/(j-i)^2

令g[i]=1/(i^2)

E[i]=∑q[j]*g[i-j]-∑q[j]*g[j-i]

∑q[j]*g[i-j]是卷積形式，求出x^i的系數就行

∑q[j]*g[j-i] (i<j<=n)

=∑q[j+i]*g[j] (i<j+i<=n)=>(0<j<=n-i)

翻轉q得到p[n-i]=q[i],t=n-i

=∑q[j+i]*g[j]=∑p[n-i-j]*g[j]=∑q[t-j]*g[j] (0<j<=t)

∑q[t-j]*g[j] (0<j<=t)是卷積形式，求出x^(n-i)的系數就行

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<algorithm>
  4 #include<cstring>
  5 #include<cmath>
  6 using namespace std;
  7 double pi=acos(-1.0);
  8 const int N=400001;
  9 struct Complex
 10 {
 11     double r,i;
 12     Complex (double r=0,double i=0):r(r),i(i){};
 
 13     Complex operator+(const Complex &x)
 14     {
 15         return Complex(r+x.r,i+x.i);
 16     }
 17     Complex operator-(const Complex &x)
 18     {
 19         return Complex(r-x.r,i-x.i);
 20     }
 21     Complex operator*(const Complex &x)
 22     {
 23         return Complex(r*x.r-i*x.i,i*x.r+x.i*r);
 24     }
 25 }A[N],B[N],C[N];
 26 int n,len;
 27 double p[N],q[N],g[N],ans[N];
 28 void rader(Complex F[],int len)
 29 {int i;
 30     int j=len/2;
 31     for (i=1;i<len-1;i++)
 32     {
 33         if (i<j)swap(F[i],F[j]);
 34          int k=len/2;
 35          while (j>=k)
 36          {
 37             j-=k;
 38             k/=2;
 39          }
 40          if (j<k) j+=k;
 41     }
 42 }
 43 void FFT(Complex F[],int len,int t)
 44 {int i,j,k;
 45     rader(F,len);
 46     for (i=2;i<=len;i=i*2)
 47     {
 48         Complex wn;
 49         wn.r=cos(t*2*pi/i);wn.i=sin(t*2*pi/i);
 50         for (j=0;j<len;j+=i)
 51         {
 52             Complex E;
 53             E.r=1;E.i=0;
 54              for (k=j;k<j+i/2;k++)
 55               {
 56                     Complex u=F[k];
 57                     Complex v=E*F[k+i/2];
 58                 F[k]=u+v;
 59                 F[k+i/2]=u-v;
 60                 E=E*wn;
 61               }
 62         }
 63     }
 64     if (t==-1)
 65      for (i=0;i<len;i++)
 66       F[i].r/=len;
 67 }
 68 int main()
 69 {int i,j;
 70     cin>>n;
 71      for (i=1;i<=n;i++)
 72      {
 73         scanf("%lf",&q[i]);
 74      }
 75      for (i=1;i<=n;i++)
 76      p[i]=q[n-i];
 77      p[0]=q[n];
 78      for (i=1;i<=n;i++)
 79      g[i]=1.0/((double)i*i);
 80      
 81      len=1;
 82      while (len<2*n) len*=2;
 83      for (i=0;i<len;i++)
 84      A[i].r=q[i],B[i].r=g[i];
 85      FFT(A,len,1);FFT(B,len,1);
 86      for (i=0;i<len;i++)
 87      C[i]=A[i]*B[i];
 88      FFT(C,len,-1);
 89      for (i=0;i<=n;i++)
 90       ans[i]=C[i].r;
 91     memset(A,0,sizeof(A));
 92     memset(B,0,sizeof(B));
 93     memset(C,0,sizeof(C));
 94     
 95       for (i=0;i<len;i++)
 96       A[i].r=p[i],B[i].r=g[i];
 97     FFT(A,len,1);FFT(B,len,1);
 98      for (i=0;i<len;i++)
 99      C[i]=A[i]*B[i];
100      FFT(C,len,-1);
101      for (i=0;i<=n;i++)
102      ans[i]-=C[n-i].r;
103     for (i=1;i<=n;i++)
104      printf("%.3lf\n",ans[i]);
105 }

[ZJOI2014]力