传送门
题外话,那个式子巨神ypy说是电磁公式之类的,反正我觉得和电磁力什么很像。
Fj=∑i<j(i−j)2qiqj−∑i>j(i−j)2qiqj
考虑化简:
第一步:
Ei=∑j<i(i−j)2qj−∑j>i(i−j)2qj
谁都会。
令f[i]=i21,g[i]=qi
有Ei=∑j<if[j]g[i−j]−∑j<if[j]g[j−i]
为什么第二个式子不能化成f[j]g[i−j],第一眼看上去,这样显然不对,那还让你做,第二眼看上去,其实有一个重要的性质不能忘记:多项式下标不能为负数,i−j<0不能作为下标。
注意到第一个式子就是卷积的形式,但是第二个式子j+(j−i)不是一个常数,考虑翻转。
令g′[i]=g[n−i],那么我们有g′[n−j+i]=g[j−i],于是式子化成Ei=∑j<if[j]g[i−j]−∑j<if[j]g′[n−j+i]
但是,做两次FFT有一点麻烦,考虑将g数组开两倍,同时包含g,g′,即对于i<n,g[i]=(n−i)21,对于i>n,g[i]=(i−n)21
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
| #include <bits/stdc++.h> #define MAXN 1000005 using namespace std; const double PI=acos(-1.0); struct Complex{ double x,y; }; inline Complex operator + (const Complex &A,const Complex &B){ return Complex{A.x+B.x,A.y+B.y}; } inline Complex operator - (const Complex &A,const Complex &B){ return Complex{A.x-B.x,A.y-B.y}; } inline Complex operator * (const Complex &A,const Complex &B){ return Complex{A.x*B.x-A.y*B.y,A.x*B.y+A.y*B.x}; } int r[MAXN]; Complex a[MAXN],b[MAXN]; inline void FFT(Complex *A,int n,int type){ for (register int i=1;i<=n;++i) if (i<r[i]) swap(A[i],A[r[i]]); for (register int i=1;i<n;i<<=1){ int R=(i<<1); Complex Wn=Complex{cos(2*PI/R),type*sin(2*PI/R)}; for (register int j=0;j<n;j+=R){ Complex w=Complex{1,0}; for (register int k=0;k<i;++k,w=w*Wn){ Complex x=A[j+k],y=A[i+j+k]*w; A[j+k]=x+y,A[i+j+k]=x-y; } } } } int main(){ int n; scanf("%d",&n); for (register int i=1;i<=n;++i){ scanf("%lf",&a[i].x); } for (register int i=0;i<n;++i){ b[i].x=-1.0/((double)(n-i)*(n-i)); } for (register int i=n+1;i<=2*n;++i){ b[i].x=1.0/((double)(i-n)*(i-n)); } int sz=2*n,m=1,L=0; while (m<=sz) m<<=1,L++; for (register int i=0;i<=m;++i){ r[i]=(r[i>>1]>>1|((i&1)<<(L-1))); } FFT(a,m,1),FFT(b,m,1); for (register int i=0;i<=m;++i) a[i]=a[i]*b[i]; FFT(a,m,-1); for (register int i=n+1;i<=2*n;++i) printf("%.3lf\n",a[i].x/m); }
|