P3338

发表于 分类于
[ZJOI2014]力

[ZJOI2014]力

差卷积的实现。

库仑力？我物理但都不会。。。

$$
\begin{aligned}
E_j=& \dfrac{F_j}{q_j}
\\
=& \sum_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n}\dfrac{q_i}{(i-j)^2}
\end{aligned}
$$

定义 $f_i=q_i$，$g_i=1/i^2$，特别地，令 $g_0=0$，于是我们有：

$$
\begin{aligned}
E_j=& \sum _ {i=0} ^ {j} f _ i g _ {j-i} - \sum _ {i=j} ^ {n} f _ i g _ {i-j}
\\
=& \sum _ {i=0} ^ {j} f _ i g _ {j-i} - \sum _ {i=0} ^ {n-j} f _ {i+j} g _ {i}
\end{aligned}
$$

令 $f ^ {\prime} _ i = f _ {n-i}$，则 $f _ {i} = f ^ {\prime} _ {n-i}$，于是我们有：

$$
E_j= \sum _ {i=0} ^ {j} f _ i g _ {j-i} - \sum _ {i=0} ^ {n-j} f ^ {\prime} _ {(n-j)-i} g _ {i}
$$

前面和后面的卷积分开做，后面那一部分的卷积做完之后需要再翻转一次然后贡献给最终数列。

时间复杂度 $O(n\log n)$。

代码：

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
55
56
57
#include<iostream>
#include<cstdio>
#include<cmath>
#define ll long long
#define ld double
using namespace std;
namespace Ehnaev{
  inline ll read() {
    ll ret=0,f=1;char ch=getchar();
    while(ch<48||ch>57) {if(ch==45) f=-f;ch=getchar();}
    while(ch>=48&&ch<=57) {ret=(ret<<3)+(ret<<1)+ch-48;ch=getchar();}
    return ret*f;
  }
}using Ehnaev::read;

const ll N=2e5;
const ld Pi=acos(-1);

struct CP{
  ld x,y;inline CP(ld a=0,ld b=0) {x=a;y=b;}
  inline CP operator+(const CP &r)const{return CP(x+r.x,y+r.y);}
  inline CP operator-(const CP &r)const{return CP(x-r.x,y-r.y);}
  inline CP operator*(const CP &r)const{return CP(x*r.x-y*r.y,x*r.y+y*r.x);}
}ans[N*2+5],f[N*2+5],g[N*2+5],h[N*2+5];

ll n,m;
ll rev[N*2+5];

inline void FFT(CP *f,bool op) {
  for(ll i=0;i<n;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
  for(ll p=2;p<=n;p<<=1) {
    ll len=p>>1;CP tG(cos(2*Pi/(ld)p),(op?-1:1)*sin(2*Pi/(ld)p));
    for(ll k=0;k<n;k+=p) {
      CP buf(1,0);
      for(ll l=k;l<k+len;l++,buf=buf*tG) {
        CP t=buf*f[l+len];f[l+len]=f[l]-t;f[l]=f[l]+t;
      }
    }
  }
  if(op) {for(ll i=0;i<n;i++) f[i].x/=(ld)n;}
}

int main() {

  n=read();
  for(ll i=1;i<=n;i++) {
    scanf("%lf",&f[i].x);h[n-i].x=f[i].x;g[i].x=(ld)(1.0/i/i);
  }
  for(m=n,n=1;n<=(m<<1);n<<=1);
  for(ll i=0;i<n;i++) rev[i]=(rev[i>>1]>>1)|((i&1)?(n>>1):0);
  FFT(f,0);FFT(g,0);FFT(h,0);
  for(ll i=0;i<n;i++) {f[i]=f[i]*g[i];h[i]=h[i]*g[i];}
  FFT(f,1);FFT(h,1);
  for(ll i=1;i<=m;i++) printf("%.3lf\n",f[i].x-h[m-i].x);

  return 0;
}