P3338 [ZJOI2014]力(FFT)

题目

P3338 [ZJOI2014]力

做法

普通卷积形式为:\(c_k=\sum\limits_{i=1}^ka_ib_{k-i}\)

其实一般我们都是用\(i=0\)开始的,但这题比较特殊,忽略那部分,其他的直接按下标存下来,反正最后的答案是不变的

好了步入正题吧,我们定义 $$F_j=\sum\limits_{i<j}\dfrac{q_iq_j}{(i-j)2}-\sum\limits_{i<j}\dfrac{q_iq_j}{(i-j)2}$$

求\(E_i=\dfrac{F_i}{q_i}\)

显然$$E_j=\sum\limits_{i<j}\dfrac{q_i}{(i-j)2}-\sum\limits_{i<j}\dfrac{q_i}{(i-j)2}$$

\[E_j=\sum_{i=1}^{j-1}\dfrac{q_i}{(i-j)^2}-\sum_{i=j+1}^{n}\dfrac{q_i}{(i-j)^2}
\]

令\(f_i=q_i,g_i=\dfrac{1}{i^2}\),特别地,\(g_0=0\),则有:$$E_j=\sum_{i=1}{j-1}f_ig_{i-j}-\sum_{i=j+1}{n}f_ig_{i-j}$$

左边部分很简单就能化成卷积形式:$$\sum_{i=1}{j-1}f_ig_{i-j}=\sum_{i=1}{j-1}f_ig_{j-i}=\sum_{i=1}^{j}f_ig_{j-i}$$

右边部分:$$\sum\limits_{i=j+1}{n}f_ig_{i-j}=\sum\limits_{i=1}{n-j}f_{i+j}g_{i}$$

令\(p_i=f_{n-i}\),则\(p_{n-i-j}=f_{i+j}\),则有:$$\sum\limits_{i=1}{n-j}f_{i+j}g_{i}=\sum_{i=1}{n-j} p_{n-i-j}g_i$$

都化成卷积形式了直接FFT

My complete code

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cmath>
using namespace std;
const double Pi=acos(-1.0);
const int maxn=300000;
struct complex{
double x,y;
complex (double xx=0,double yy=0){
x=xx,y=yy;
}
}f[maxn],p[maxn],g[maxn];
int n,limit,L;
int r[maxn];
complex operator + (complex x,complex y){
return complex(x.x+y.x,x.y+y.y);
}
complex operator - (complex x,complex y){
return complex(x.x-y.x,x.y-y.y);
}
complex operator * (complex x,complex y){
return complex(x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x);
}
inline void FFT(complex *A,int type){
for(int i=0;i<limit;++i)
if(i<r[i])
swap(A[i],A[r[i]]);
for(int mid=1;mid<limit;mid<<=1){
complex WN(cos(Pi/mid),type*sin(Pi/mid));
for(int R=mid<<1,j=0;j<limit;j+=R){
complex w(1,0);
for(int k=0;k<mid;++k,w=w*WN){
complex x=A[j+k],y=w*A[j+mid+k];
A[j+k]=x+y,
A[j+mid+k]=x-y;
}
}
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;++i){
scanf("%lf",&f[i].x),f[i].y=0,
g[i].x=(double)1.0/i/i;
p[n-i].x=f[i].x;
}
limit=1,L=0;
while(limit<n+n)
limit<<=1,
++L;
for(int i=0;i<limit;++i)
r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
FFT(f,1),FFT(p,1),FFT(g,1);
for(int i=0;i<limit;++i)
f[i]=f[i]*g[i];
for(int i=0;i<limit;++i)
p[i]=p[i]*g[i];
FFT(f,-1),FFT(p,-1);
for(int i=1;i<=n;++i)
printf("%.3lf\n",f[i].x/limit-p[n-i].x/limit);
return 0;
}/*
*/
上一篇:给jdk写注释系列之jdk1.6容器(7)-TreeMap源码解析


下一篇:深入理解CADisplayLink和NSTimer