Steven_MengのBlog

传送门

题外话,那个式子巨神ypy说是电磁公式之类的,反正我觉得和电磁力什么很像。

$F_j=\sum_{ij}\dfrac{q_i q_j}{(i-j)^2}$

考虑化简:

第一步:

$E_i=\sum_{ji}\dfrac{q_j}{(i-j)^2}$

谁都会。

令$f[i]=\dfrac1}{i^2},g[i]=q_i$

有$E_i=\sum _{j<i} f[j]g[i-j]-\sum _{j<i} f[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]$,于是式子化成$E_i=\sum _{j<i} f[j]g[i-j]-\sum _{j<i} f[j]g’[n-j+i]$

但是,做两次FFT有一点麻烦,考虑将$g$数组开两倍,同时包含$g,g’$,即对于$in,g[i]=\frac{1}{(i-n)^2}$\dfrac

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);
}

 评论